PaperlayerView exact revision on arXiv

Source-linked AI summary

Enzo+Moray: Radiation Hydrodynamics Adaptive Mesh Refinement Simulations with Adaptive Ray Tracing

John H. Wise, Tom Abel

arXiv:1012.2865v2astro-ph.IMastro-ph.CO

TL;DR

Multidimensional radiation-hydrodynamics simulations require a method that captures radiation transport, gas response, and changing spatial resolution. The paper develops photon-conserving adaptive ray tracing coupled to Enzo’s AMR hydrodynamics, chemistry, and energy solvers, with tests showing accurate benchmark behavior and scalable parallel execution. The method is intended for dynamical astrophysical problems including stellar and black-hole feedback, while retaining limitations from timestep size, sampling, and approximation transitions.

  • Problem

    Multidimensional radiation calculations are difficult because radiation depends on seven variables and has nonlocal thermal and hydrodynamical effects.

  • Method

    The paper couples photon-conserving adaptive ray tracing with Enzo’s AMR hydrodynamics, chemistry, and energy solvers, using geometric corrections and an optically thin multiple-source approximation.

  • Results

    The implementation reproduces expected analytical behavior with three photon energy groups for Population III stellar radiation and reaches a late-time ionisation-front radius 4% above the Strömgren radius in a benchmark test.

  • Takeaways & Limitations

    Enzo+Moray provides a radiation-hydrodynamics framework for studying dynamical stellar, black-hole, and magnetised-gas problems within an AMR cosmological code.

  • Takeaways & Limitations

    Large timesteps can slow ionisation-front propagation, while the optically thin approximation can produce angular spike artifacts near its transition to full ray tracing.

Abstract

from arXiv · show

We describe a photon-conserving radiative transfer algorithm, using a spatially-adaptive ray tracing scheme, and its parallel implementation into the adaptive mesh refinement (AMR) cosmological hydrodynamics code, Enzo. By coupling the solver with the energy equation and non-equilibrium chemistry network, our radiation hydrodynamics framework can be utilised to study a broad range of astrophysical problems, such as stellar and black hole (BH) feedback. Inaccuracies can arise from large timesteps and poor sampling, therefore we devised an adaptive time-stepping scheme and a fast approximation of the optically-thin radiation field with multiple sources. We test the method with several radiative transfer and radiation hydrodynamics tests that are given in Iliev et al. (2006, 2009). We further test our method with more dynamical situations, for example, the propagation of an ionisation front through a Rayleigh-Taylor instability, time-varying luminosities, and collimated radiation. The test suite also includes an expanding H II region in a magnetised medium, utilising the newly implemented magnetohydrodynamics module in Enzo. This method linearly scales with the number of point sources and number of grid cells. Our implementation is scalable to 512 processors on distributed memory machines and can include radiation pressure and secondary ionisations from X-ray radiation. It is included in the newest public release of Enzo.

1 INTRODUCTION

Radiative transfer is important across astrophysical systems but difficult to model in multiple dimensions because radiation depends on seven variables and responds nonlocally. The paper introduces Enzo+Moray to couple adaptive ray tracing with AMR cosmological hydrodynamics, extending simulations beyond static-density treatments.

  • Multidimensional radiation calculations are difficult because the radiation field depends on three spatial, two angular, frequency, and time variables, while thermal and hydrodynamical responses are nonlocal.
  • Radiation heats and ionises gas, producing expanding H II regions whose fronts can transition from supersonic R-type to subsonic D-type propagation.
  • Moment methods are fast and independent of source number, whereas ray-tracing methods propagate radiation along rays through computational grids or particle sets.
  • Earlier three-dimensional calculations often used static density fields and therefore omitted the hydrodynamical response to radiation.
  • Enzo+Moray implements adaptive ray tracing in Enzo’s cosmological hydrodynamics AMR framework and couples the radiation field to hydrodynamics at small timescales.
  • The implementation includes geometric corrections, additional radiation physics, and an optically thin multiple-source approximation to improve radiative-transfer calculations.

2 TREATMENT OF RADIATIVE TRANSFER

The method solves radiative transfer from point sources with adaptive ray tracing, integrating photon absorption along ray segments and coupling the radiation field to numerical calculations. It improves angular sampling with ray splitting and geometric corrections while retaining time-dependent propagation and accounting for known conservation trade-offs.

  • Adaptive Ray Tracing: The required ray count grows as N_ray = 4πR^2/(Δx)^2, making uniform sampling at large radii computationally wasteful.At smaller radii, the same rays oversample cells relative to the target radius, motivating adaptive splitting.
  • Adaptive Ray Tracing: Adaptive ray tracing progressively splits HEALPix-based rays when angular sampling becomes too coarse.Initial rays distribute each source’s photon luminosity equally across HEALPix pixels; child rays inherit directionally refined sampling and one quarter of the parent photon flux.
  • Ray Propagation: Each ray is propagated until its radiative-transfer timestep ends, its flux is more than 99.9% absorbed in one cell, or it exits the domain.Rays halted at a timestep are saved for the next step, preserving changing hydrodynamical states, ionisation states, and luminosities.
  • Radiation Field: Photon-conserving rates are obtained by integrating the transfer equation along ray segments and equating photo-ionisation to absorption.Photo-ionisation and photo-heating contributions from individual rays are summed in each cell.
  • Geometric Corrections: A geometric overlap correction suppresses non-spherical artifacts caused when a ray’s angular pixel does not cover an entire cell face.The correction multiplies photo-ionisation and photo-heating rates by the overlap factor fc, while absorbed radiation remains unchanged.
  • Geometric Corrections: The geometric correction sacrifices exact photon conservation in its implementation, although the authors report no perceptible deviations in tested cases.A fully conservative treatment would also account for neighboring cells overlapped by the ray’s solid angle.
  • Secondary Ionisations: In gas with low electron fractions, most deposited photon energy produces hydrogen and helium ionisations, whereas nearly ionised gas receives mostly photo-heating.The energy-partition factor approaches one as the electron fraction approaches one.

3 NUMERICAL IMPLEMENTATION IN ENZO

The implementation integrates adaptive ray tracing into Enzo’s parallel block-structured AMR framework and includes radiation-hydrodynamic coupling, timestep control, and distributed-memory ray tracing. It was included in Enzo v2.0.

  • Implementation: The implementation adds adaptive ray tracing to Enzo, a publicly available parallel block-structured AMR code.The section covers photon-package handling, radiation hydrodynamics, timestep selection, and ray tracing across the AMR hierarchy.
  • Radiation Hydrodynamics: Radiation is coupled to the hydrodynamics solver at small timescales for radiation-hydrodynamic calculations.
  • Release: The implementation was included in the v2.0 public version of Enzo.

3.1 Programming Design

The module creates and traces photon packages across AMR grids, adaptively splitting rays and coupling radiative updates to chemistry, energy, and hydrodynamics.

  • Algorithm overview: Photon packages are created from point sources on an initial HEALPix level and placed in the highest-resolution AMR grid containing each source.Each package stores flux, type, energy, timing, radius, column density, HEALPix information, and source position.
  • Ray traversal: The solver loops over AMR grids and traces each ray by computing its direction, position, cell crossing, and segment length.Rays are moved between grids when subgrids or grid boundaries are encountered, with MPI communication for remote destination grids.
  • Design constraint: The implementation cannot use stored ray-segment paths because grid structure and point-source positions are unrestricted and may change.This removes an available computational acceleration strategy.
  • Adaptive ray tracing: Rays split into four child rays when their solid angle exceeds the user-defined sampling criterion.The parent ray is deleted after its children are inserted into the linked list.
  • Radiative updates: Each ray updates geometric correction, optical depth, photo-ionisation, photo-heating, column density, time, flux, and radius while traversing cells.Optional modules add secondary X-ray ionisations, Compton heating, and radiation pressure; rays terminate when exhausted or optically thick.
  • Solver coupling: After radiation fields are updated, the chemistry and energy solver advances affected cells before the photon time is advanced through the AMR hierarchy.The overall flow includes photon creation, ray tracing, transport between AMR grids, and hydrodynamic coupling.

3.2 Energy groups

Enzo+Moray represents radiation with monochromatic photon packages grouped by photon type, trading additional rays for early termination of strongly absorbed components.

  • Energy-group representation: Photon packages are monochromatic energy groups assigned to hydrogen, helium, X-ray, or molecular-hydrogen-dissociating radiation types.The implementation distinguishes H I, He I, He II, X-ray, and Lyman-Werner photons.
  • Trade-off: The monochromatic approach increases ray counts with the number of frequency bins but permits early termination of fully absorbed rays.This can reduce work for components with high absorption cross-sections or low initial intensity.
  • Energy-group choice: One energy group per photon type is sufficient to match expected analytical tests.For Population III stellar radiation, the implementation uses three groups—H I, He I, and He II—with energies equal to average photon energies above threshold.

3.3 Coupling with Hydrodynamics

Radiation is coupled to Enzo through frequent radiation updates, targeted chemo-thermal subcycling, and the existing operator-split gravity, hydrodynamics, and chemistry sequence.

  • Motivation: Radiative transfer must update the gas frequently because radiation can change on shorter timescales than ordinary hydrodynamical evolution.Without intermediate updates, an ionisation front can appear to enter a grid only at the next hydrodynamical timestep.
  • Operator-split sequence: On the finest AMR level, Enzo solves radiation first, then updates species fractions and energies for irradiated cells with a non-equilibrium chemistry solver.The main grid loop subsequently handles gravity, hydrodynamics, non-irradiated chemistry, particle motion, and star formation.
  • Chemistry and energy solver: The chemistry and energy solver uses subcycles limited to a 10% change in electron density, neutral hydrogen density, or specific energy.The subcycle timestep also depends on the hydrodynamic timestep and electron density.
  • Selective updates: The implementation dynamically separates cells with and without radiation, updating irradiated cells on radiation timesteps and the remainder on hydrodynamic subcycles.This targets expensive chemo-thermal updates to cells affected by radiation.
  • Time-state mismatch: When cells first acquire radiation, their initial chemo-thermal state can lag the radiation transport time, but the authors omit a correction because relevant unirradiated-gas timescales are short.They report no resulting inaccuracies in the cited test problems and applications.

3.4 Temporal evolution

Temporal control balances accuracy against computational cost by limiting neutral-fraction changes, smoothing noisy timestep estimates, and using physically motivated constraints when ionisation-front speed matters.

  • Timestep accuracy: Large radiation timesteps make ionisation fronts advance too slowly because the chemo-thermal state remains fixed while rays traverse cells.For a single H II region, the front speed is limited to approximately λ_mfp/dtP.
  • Neutral-fraction control: The adaptive scheme restricts neutral-fraction changes within the ionisation front and ignores ionised cells whose absolute changes are small.The implementation considers cells with optical depth τ > 0.5 when determining the timestep.
  • Timestep stabilization: Spatial Gaussian smoothing and a running average stabilize the minimum cell-based timestep, whose raw values are noisy across space and time.Figure 4 compares the noisier unmodified timestep with the smoother running-average result while limiting neutral-fraction change to 5%.
  • Step-size limiter: The timestep is constrained to less than twice the previous timestep to prevent inaccurate jumps caused by overly large steps.The resulting evolution of dtP is smoother than the unsmoothed minimum cell values in the growing Strömgren sphere.
  • Alternative strategy: Time-averaged quantities can preserve analytical agreement with 10^-3 fewer timesteps than a method without time averaging, but this approach was still being implemented.The method iteratively averages neutral density, optical depth, photo-ionisation rate, and electron density until convergence.
  • Physically motivated timesteps: A large timestep may be acceptable during a D-type phase but can make the ionisation front lag when it later accelerates and becomes R-type.The two stages can differ in velocity by up to a factor of ~10, while the numerical solution eventually approaches the correct radius.

3.5 Parallelisation Strategy

The parallelisation strategy distributes photon packages across processors with non-blocking MPI communication, while testing shows ray-level OpenMP threading scales better than partitioning work over grids.

  • MPI communication: Photon packages are grouped by destination processor and exchanged with non-blocking MPI to overlap communication and computation.Messages are buffered, posted asynchronously, and processed as they arrive.
  • MPI communication: Each MPI message carries at most Nmax = 10^5 rays, with receive buffers allocated accordingly.The implementation counts destination messages and communicates their number before transferring photon data.
  • Ray transport: Rays moved to local grids are transported immediately, while remote rays are received through non-blocking MPI calls.The algorithm checks for already-arrived messages and waits only when necessary.
  • Termination: Ray transport completes when all processors exhaust their workloads and rays are absorbed, leave the domain, or travel cdtP.Completion is checked using a non-blocking procedure analogous to the message-count exchange.
  • Hybrid parallelisation: OpenMP/MPI experiments found that partitioning work over grids scales poorly, whereas threading over rays within each grid is more effective.Because rays use linked lists, the implementation manually splits each list among threads.

4 RADIATIVE TRANSFER TESTS

The radiative-transfer tests evaluate ionisation-front expansion, temperature evolution, dense-clump shadowing, and multiple-source cosmological reionisation. Across these cases, the method reproduces analytical expectations and comparison-code results while exposing controlled ray-splitting artifacts.

  • Test design: The test suite follows the Cosmological Radiative Transfer Codes Comparison Project and includes problems designed to reveal asymmetries, artifacts, and implementation errors.The tests cover standard radiative-transfer cases and more complex hydrodynamical configurations.
  • Test 1: Isothermal expansion: In Test 1, the ionisation-front width is ∼0.7 kpc, matching the expected 18λmfp = 0.74 kpc for 13.6 eV monochromatic radiation.The calculated front radius agrees closely with the analytical Strömgren solution.
  • Test 1: Isothermal expansion: Test 1 preserves spherical symmetry and agrees with the more accurate analytical solution, while ray-splitting artifacts produce small radial plateaus.The artifacts occur near 1.5 and 3 kpc and do not affect the overall solution.
  • Test 2: Temperature evolution: For the harder blackbody spectrum, the ionisation front is thicker and the total neutral fraction decreases to 0.67 over 4trec = 500 Myr.The result agrees with the analytical expectation and other RT06 codes.
  • Test 2: Temperature evolution: In the blackbody test, rIF initially lags Rs by 10%, then asymptotes to a solution 4% greater than Rs.Using one energy group with mean energy 29.6 eV gives rIF/Rs = 1.08.
  • Test 3: Dense-clump shadowing: The dense-clump test produces a partially trapped ionisation front, gradual increases in ionised fraction and temperature, and sharp shadows characteristic of non-diffusive ray tracing.A discretisation artifact remains because the clump surface was not smoothed.
  • Test 4: Multiple sources: In the multiple-source cosmological test, the volume-averaged ionised fraction overtakes the mass-averaged fraction at t ∼170 kyr, and 65% of the simulation is ionised at the end.The ordering indicates inside-out reionisation, while photo-heated regions extend 2–3 times farther than ionised regions.

5 RADIATION HYDRODYNAMICS TESTS

The radiation-hydrodynamics tests reproduce established ionisation-front behavior while exposing how front definitions, density structure, and gas dynamics shape the results. They also demonstrate photo-evaporation and symmetry effects in dense-clump problems.

  • Test suite: The suite tests H II-region expansion in uniform media, isothermal spheres, and dense clumps against analytical or published radiation-hydrodynamics results.The uniform-medium tests use analytical comparisons, while the isothermal-sphere test is compared with RT09 because no analytical full-gas-dynamics solution exists.
  • Uniform H II region: The uniform-medium calculation does not capture pressure equilibrium at rf = 185 kpc, but it does capture the R-type-to-D-type transition near Rs = 5.4 kpc.The test domain is therefore sufficient for the transition but not the eventual equilibrium radius.
  • Uniform H II region: Before 2trec ≈250 Myr, the T = 10^4 K front definition overestimates the analytical radius by over 10%, whereas xe = 0.5 underestimates it by ∼20%.The temperature criterion matches the late-time t^4/7 growth more closely, while the ionised-fraction criterion remains systematically smaller.
  • Uniform H II region: The uniform-medium H II region drives a forward shock, while high-energy photons partially ionise and heat exterior gas, producing an inward photo-evaporative flow.The shocked gas can recombine without radiatively cooling, contributing to the broad ionisation-front structure.
  • Isothermal sphere: The isothermal-sphere front agrees with RT09 within 5%, transitioning from a weak R-type front to a shock-producing phase as it moves through the density gradient.The inner density decreases by over two orders of magnitude after 25 Myr; reflecting boundaries create a slight shock-broadening artifact.
  • Dense clump: Photo-evaporation of a dense clump produces an over-pressurised heated layer that expands into the ambient medium and forms photo-evaporative flows.By 50 Myr, only the core and its shadow remain neutral, while the core is compressed by the surrounding warm medium.

6 RADIATION HYDRODYNAMICS APPLICATIONS

The applications extend the framework to asymmetric, unstable, time-dependent, and magnetised radiation-hydrodynamics scenarios. They show that radiation interactions generate champagne flows, instability-driven structure, source-history signatures, and magnetic confinement.

  • Scope: The applications cover champagne flows, irradiated Rayleigh-Taylor instability, irradiated blast waves, collimated radiation, variable luminosity, and expanding H II regions in magnetic fields.These cases extend the standard RT06 and RT09 test suites to more dynamical configurations.
  • Champagne flow: After escaping an overdense clump, the champagne-flow ionisation front becomes parabolic and perturbations grow through Rayleigh-Taylor instability during accelerated outflow.The perturbations originate from mismatches between HEALPix pixels and the Cartesian grid, then grow as gas exits the clump.
  • Rayleigh-Taylor instability: In the optically thick Rayleigh-Taylor case, a D-type front, shock reflections, and photo-evaporation transform the initial structure into a turbulent medium.The instability develops a chaotic jet-like structure that is eventually destroyed by photo-evaporation and interactions with the hot medium.
  • Rayleigh-Taylor instability: The optically thin Rayleigh-Taylor case is less violent because its R-type front interacts less strongly with the initial instability.Heating produces a warm, partially ionised turbulent medium while the dense gas remains optically thick.
  • Variable luminosity: A variable source luminosity leaves the overall H II-region growth nearly unchanged while producing transient ∼500 pc density perturbations that smooth over ∼50 Myr.The photo-ionisation-rate shells retain the source’s luminosity history and decline as 1/r^2 when averaged over equal-width shells.
  • MHD application: In a magnetised medium, the expanding H II region becomes oblate after t > 0.5 Myr because magnetic pressure confines expansion perpendicular to the field.The gas can then flow primarily along magnetic field lines.

7 RESOLUTION TESTS

Resolution tests show that photon conservation preserves accurate front propagation under marginal spatial resolution, while angular, frequency, and timestep choices control radiation-field smoothness and reactive-flow details.

  • Spatial resolution: At 16^3 and 32^3 resolution, the ionisation front is initially underestimated but converges within 0.5% of higher-resolution runs.The discrepancy persists for roughly the first 25 Myr before convergence.
  • Spatial resolution: The 16^3 and 32^3 runs lag by 1% until 300 Myr and are later larger by 0.5% than higher-resolution cases.The authors attribute this accuracy in marginally resolved cases to photon conservation.
  • Angular resolution: For angular sampling Φc ≤ 3.1, insufficient rays create cell-to-cell photo-ionisation-rate variations, while higher Φc progressively removes shell artifacts.At Φc = 5.1, cell variations disappear but a shell artifact remains; higher sampling removes it as well.
  • Angular resolution: The deviation from an inverse-square radiation field decreases as σ ∝ Φc^-0.6, independently of spatial resolution.Angular resolution is controlled relative to cell widths rather than absolute solid angles.
  • Frequency resolution: The ionisation-front solution converges for nν ≥ 4, while monochromatic radiation misses the double-peaked density and velocity structure produced by multi-frequency spectra.Higher-energy photons penetrate to larger column densities and heat gas farther from the source.
  • Temporal resolution: Variable timesteps based on dnH/dt or dI/dt agree within a few percent throughout the simulation, whereas larger constant timesteps lag behind the analytical solution.The adaptive methods shorten timesteps during rapid early expansion and lengthen them as the front approaches the Strömgren radius.

8 METHODOLOGY TESTS

The methodology tests evaluate geometric corrections, optically thin approximations, high-energy secondary ionisations, radiation pressure, and multiple-source radiation fields. Results show improved angular symmetry and physically consequential radiation effects, while the optically thin transition remains artifact-prone.

  • 8.1 Geometric correction: Geometric correction reduces shell artifacts and improves spherical symmetry, although a 45-degree non-spherical artifact is only partially removed.The correction accounts for partial solid-angle coverage of cells, adding geometric dilution beyond ray splitting.
  • 8.2 Optically-thin approximation: The optically thin approximation produces angular spikes at the τthin = 0.1 interface because cells can receive inconsistent thin and thick ray treatments.The authors had not determined a technique to avoid these artifacts.
  • 8.3 X-Ray secondary ionisations and reduced photo-heating: Secondary ionisations from a 1 keV spectrum expand the hot T = 10^5 K core by a factor of 2 and increase outer-region ionisation by ∼10, lowering temperatures by a factor of 2–3.High-energy photons distribute more energy into hydrogen and helium ionisations, leaving less energy for photo-heating.
  • 8.4 Radiation pressure: Radiation pressure increases the ionisation-front radius by ∼5% = 0.16 pc at 140 kyr and raises maximum ionised-gas velocity from 10 km s−1 to 50 km s−1.It also reduces central density by a factor of 20 at 40 kyr and gas density by a factor of 10 at 140 kyr.
  • 8.5 Consolidated H ii region with two sources: With two sources, individual H ii regions remain spherically symmetric before overlap and form a common elongated region without perceptible multiple-source artifacts.The regions grow to r = 4 kpc before overlapping.

9 PARALLEL PERFORMANCE

Parallel tests assess weak and strong scaling of enzo+moray. The code scales well in many components and overall to 512 processes, but ray communication and load imbalance limit performance in some configurations.

  • 9.1 Weak Scaling: Weak scaling uses one 64^3 block per process with a source per block and measures performance at constant work per processor.The test uses a uniform, non-AMR grid and runs for 250 Myr.
  • 9.1 Weak Scaling: Above 8 processes, most code components show good weak scaling, while inter-processor ray communication degrades and becomes dominant above 512 processes.The ray communication routines exhibit an N^1.5 dependence, whereas other components scale well.
  • 9.2 Strong Scaling: Strong scaling uses a 256^3 AMR cosmological-reionisation calculation with ∼392^3 zones, where hydrodynamics and non-equilibrium chemistry scale well.Load balancing becomes problematic at 128 and 256 cores, especially for grids containing multiple radiation sources.
  • 9.2 Strong Scaling: Ray communication does not scale in the strong-scaling problem because grids with multiple sources create unequal ray-tracing workloads and processor idle time.Splitting AMR grids by ray-tracing work could help, but may increase boundary-related overhead elsewhere in Enzo.

10 SUMMARY

The paper presents enzo+moray as a photon-conserving radiation-hydrodynamics implementation coupling adaptive ray tracing to Enzo’s AMR hydrodynamics. Tests show accurate low-resolution solutions and broad physical flexibility, while scaling reaches O(10^3) processors for large reionisation problems.

  • 10 SUMMARY: enzo+moray couples photon-conserving adaptive ray tracing to Enzo’s cosmological AMR hydrodynamics, forming a fully functional radiation hydrodynamics code.The method can produce accurate solutions with coarse spatial resolution.
  • 10 SUMMARY: A geometric correction factor improves ray tracing on Cartesian grids, and the code remains accurate at low spatial, angular, and frequency resolution.The stated exception is the use of large constant timesteps.
  • 10 SUMMARY: The implementation supports analytical, radiative-transfer comparison, and dynamical tests including champagne flows, Rayleigh-Taylor instabilities, beamed radiation, time-varying sources, and MHD H ii regions.These tests demonstrate the flexibility and fidelity of enzo+moray across varied radiation-hydrodynamics problems.
  • 10 SUMMARY: The code scales well to O(10^3) processors for problems with ∼10^9 computational cells and ∼10^4 sources, while also showing good strong scaling for sufficiently large AMR calculations.The full implementation is included in the public Enzo release.

Neutral Fraction

The supplied passages contain repeated numeric labels and a “Ray-to-cell sampling ( c)” label, but no interpretable neutral-fraction result.

  • The sequence “0 5 10 15 20” is repeated across four passages.
  • No passage states a neutral-fraction value, trend, comparison, or conclusion.
  • “Ray-to-cell sampling ( c)” is the only distinct textual label supplied.
Loading 1012.2865v2…