Enabling Radiative Transfer on AMR grids in CRASH

12.9 MB·

by  N. Hariharan

MNRAS000,1–17(2015) Preprint18January2017 CompiledusingMNRASLATEXstylefilev3.0
Enabling Radiative Transfer on AMR grids in CRASH
N. Hariharan1,4(cid:63), L. Graziani1,5, B. Ciardi1, F. Miniati2, H.-J. Bungartz3
1Max-Planck-Institu¨t fu¨r Astrophysik, Karl-Schwarzschild-Strasse 1, 85748 Garching, Germany
2Department of Physics, ETH-Zu¨rich, Wolfgang-Pauli-Strasse 27, CH-8093, Zu¨rich, Switzerland
3Institu¨t fu¨r Informatik, TU Mu¨nchen, D-80290, Mu¨nchen, Germany
4Intel Technology India Private Limited, #136 Old Airport Road, Bangalore-560008, India
5INAF Osservatorio Astronomico di Roma, Via Frascati 33, 00040, Monte Porzio Catone (RM), Italy
7
1 Accepted1988December15.Received1988December14;inoriginalform1988October11
0
2
Key words: Cosmology:theory—Radiative Transfer—methods:numerical
n
a
J

6
1 ABSTRACT

]
We introduce CRASH-AMR, a new version of the cosmological Radiative Transfer
O
(RT) code CRASH, enabled to use refined grids. This new feature allows us to attain
C
higher resolution in our RT simulations and thus to describe more accurately ionisa-
. tion and temperature patterns in high density regions. We have tested CRASH-AMR by
h
p simulatingtheevolutionofanionisedregionproducedbyasinglesourceembeddedin
- gasatconstantdensity,aswellasbyamorerealisticconfigurationofmultiplesources
o in an inhomogeneous density field. While we find an excellent agreement with the
r
t previousversion of CRASH whenthe AMR featureisdisabled, showingthatnonumer-
s ical artifact has been introduced in CRASH-AMR, when additional refinement levels are
a
used the code can simulate more accurately the physics of ionised gas in high density
[
regions. This result has been attained at no computational loss, as RT simulations on

1 AMR grids with maximum resolution equivalent to that of a uniform cartesian grid
v can be run with a gain of up to 60% in computational time.
9
2 1 INTRODUCTION 2006), while helium reionisation is believed to be complete
4 at z ∼ 2.7 (e.g. Madau & Meiksin 1994; Compostella et al.
4 The observed universe shows a large variation in structures 2013).
0 aswemovealongdifferentscales.Independentobservations
A large number of simulations have been done to un-
. of the distribution of faint radio sources, optically selected
1 derstand the formation of the first structures and the sub-
galaxies and the X-ray background show that our universe
0 sequent galaxy formation and evolution. Examples of colli-
is homogeneous on a scale larger than 200 Mpc (see e.g.
7 sionlessdark-mattersimulationsincludeShawetal.(2006),
1 Wu et al. 1999 for a review), while at smaller scales it ap- Boylan-Kolchinetal.(2009),Kimetal.(2011),Anguloetal.
: pears inhomogeneous due to the presence of a large num-
v (2012)andHarnois-D´erapsetal.(2013),whilee.g.Springel
ber of structures: for example, galaxy clusters and galax-
i (2005), Vogelsberger et al. (2014) and Schaye et al. (2014)
X ies at tens of Mpc, or stellar clusters and molecular clouds
follow also the gas dynamics. To study radiative feedback,
at parsec scales. In the ΛCDM concordance model of our
r radiativetransfer(RT)simulationsareusuallyusedaspost-
a universe, the presence of primordial density perturbations
processing tools for hydro codes or N-body codes, and they
led to the gravitational collapse and cooling of gas in pre-
form a major aspect of the study of structure formation in
existingdark-matterhalos,andtothesubsequentformation
general, and the IGM reionisation in particular (e.g. Ciardi
of radiating sources like stars and quasars (e.g. Mo et al.
etal.2003;Ilievetal.2006a,2007;Zahnetal.2007;Trac&
2010). The chemical, mechanical and radiative feedback by
Cen 2008; Ahn et al. 2012; Iliev et al. 2014; Graziani et al.
thesesourcesontheirsurroundingsinducesacomplexinter-
2015). Simulations where the gas dynamics and RT effects
playbetweenthegalaxyformationprocessandtheevolution
areself-consistentlyaccountedforhavealsobeenperformed
of the intergalactic medium (IGM; Ciardi & Ferrara 2005;
(e.g. Gnedin & Ostriker 1997; Semelin et al. 2007; Gnedin
Dav´e 2005; Barkana & Loeb 2007; Meiksin 2009). Among
2014; Gnedin & Kaurov 2014; Pawlik et al. 2015). Large
thesefeedbackeffects,aprominentplaceisoccupiedbythe
scale hydro and RT simulations require a large amount of
IGMreionisationprocess,whichdenotesthetransitionfrom
computational and memory resources to handle spatial res-
aneutralintergalacticgastoanIGMwhichis(almost)fully
olutionsspanningseveralordersofmagnitude.Inthecaseof
ionised in its hydrogen component by z ∼ 6 (e.g. Fan et al.
grid-basedcodes,thisproblemhasbeenaddressedtoagreat
(cid:13)c 2015TheAuthors
2 Hariharan N. et al.
extent by the use of Adaptive Mesh Refinement (AMR) point sources. A final example is IFT (Alvarez et al. 2006),
schemes (e.g. Kravtsov et al. 1997; Teyssier 2002; Miniati which has been developed to explicitly follow the I-front
& Colella 2007; Bryan et al. 2014), which provide a great around a point source.
amount of theoretical and algorithmic resources, and have Fromtheaboveexamples,itisclearthatAMR-enabled
been used extensively to solve a wide range of problems, stand-aloneRTcodesareincreasinglyinuseandtheystand
such as numerical relativity, global weather and nuclear fu- to benefit from the advantages that AMR provides. Keep-
sion modelling (e.g. Plewa et al. 2003). ing this in mind, we have developed CRASH-AMR, a novel
AMRschemesuseameshoragridtodescribethephys- implementation of the RT code CRASH (Ciardi et al. 2001;
icaldomainandtoprogressivelyincreasethegridresolution Maselli et al. 2003; Maselli & Ferrara 2005; Maselli et al.
inconfinedregionsofthemesh,basedonasetofrefinement 2009; Pierleoni et al. 2009; Partl et al. 2011; Graziani et al.
criteria. By selectively increasing the resolution only in the 2013),whichisinterfacedwiththeopensourceAMRlibrary
interesting part of the domain, AMR methods optimise the CHOMBO (Adams et al. 2011) to perform RT simulations on
global memory and computational resource requirements. AMR grids. The aim of this paper is to introduce this new
The original idea of AMR was to introduce a finer grid in version, to discuss its implementation details and to show
regionsofhighernumericalerror,whichcanbeidentifiedby theteststhatwehavedonetoverifyandvalidateCRASH-AMR.
densitygradientsorbyRichardsonextrapolation(Berger& The paper is structured as follows. In Section 2 we
Oliger 1984; Berger & Colella 1989). However, AMR can brieflyintroduceAMRanditsdifferentschemes,theCHOMBO
be employed very flexibly, and different refinement criteria libraryisalsopresentedinthiscontextandwementionsome
can be chosen. For example, in cosmological applications it of its applications. In Section 3 we introduce the ray trac-
is customary to adopt a Lagrangian criterion in which cells ing implementation in the CRASH code. In Section 4 we dis-
are refined when their mass is above a given threshold (e.g. cussthemethodusedtocoupleCRASHandCHOMBOtoenable
Wise&Abel2011).Alternatively,itisalsopossibletorefine AMR in CRASH. The tests performed to verify our code are
specific volumes in a Eulerian style, which is advantageous discussedinSection5andtheresultsaresummarisedinSec-
for example if one is interested in resolving turbulent mo- tion6.Theadvantagesofthecodeintermsofcomputational
tions inside cosmological halos (e.g. Miniati 2014). costs is highlighted in Appendix A and B. Here, we discuss
AMR is also very flexible from an algorithmic point the tests done to show the dependence on grid-resolution
of view. For example, it naturally supports divergence- and performance of the code. We discuss code performance
free magnetohydrodynamics (e.g. Miniati & Martin 2011; in terms of run-time and correctness of results when com-
Teyssier et al. 2006; Lee & Deane 2009), it has been used pared to running RT simulations on uniform grids.
to make detailed studies of thermonuclear flashes in FLASH
(Fryxell et al. 2000), and also to simulate the formation of
large-scale cosmological structures in ENZO (see e.g. Bryan 2 BASICS OF AMR AND THE CHOMBO
et al. 2014 and Plewa et al. 2003 for more examples). LIBRARY
Many hydro codes make use of AMR schemes to carry
out gas dynamic simulations over a range of spatial scales. AMR is a successful technique when a problem presents
Some of them are also coupled with RT schemes to per- a highly inhomogeneous spatial distribution implying that
form self-consistent simulations where the RT feedback someregionsrequireadditionalresolution,i.e.additionalre-
is accounted for in the dynamical evolution of the gas. finement. An AMR mesh is defined as structured (SAMR)
TheHydrodynamicsAdaptiveRefinementTree(HART)code if its cells (typically cartesian) are connected in a regu-
(Gnedinetal.2009)usestheOTVETapproximationforthe lar geometry. Un-structured grids composed of triangular
3DRTimplementation(Gnedin&Abel2001),whilethecos- or tetrahedral cells without a regular connectivity can also
mologicalhydrodynamicscodeRAMSES(Rosdahletal.2013), bedefined(seeforexampleMavriplis1997,Khokhlov1998,
designedforsimulationsofstructureformation,incorporates Springel 2011 and Paardekooper et al. 2008). We limit our
RTusingtheM1closureformalism(Levermore1984).Both discussion to SAMR schemes in this paper.
HART and RAMSES make use of the the Fully Threaded Tree Therearemanyalternativewaystorefineacertainsub-
described in Khokhlov (1998) to implement AMR. For the domain using different SAMR schemes. In the Cell-based
RT ENZO uses an adaptive ray-tracing scheme implemented (CBAMR)schemeeachcellofthegridisrefinedasandwhen
in the HEALPix library (Abel & Wandelt 2002; Go´rski et al. requiredandgenerallyaquad-tree(2D)orandoct-tree(3D)
2005).TheinterestedreadercanfindmoreexamplesinIliev formsthehierarchicalstructurerelatingthe“coarser”(par-
et al. (2009). ent)and“finer”(child)cells(Youngetal.1991).TheBlock-
EnablingAMRinastand-aloneRTcodemakesitsuit- structured AMR(BSAMR)schemetagsandrefinesaregion
able for post-processing the output of many grid-based hy- of the grid, by some integer factor, based on pre-assigned
dro codes that use the same AMR logic. By representing criteria (Berger & Oliger 1984). The various refinements
theregionsofinterestwithhighresolutiongrids,wecanfol- aretypicallyorganisedinahierarchicalstructureconnecting
low the details of the growth of the ionised bubbles around coarserlayers(parentlevel)withrefinedlayers(childlevel);
the sources and understand the impact of the RT feedback the grids are then stored and maintained independent from
effectsonthesurroundingenvironment.Anumberofstand- eachother(Berger&Oliger1984).Whenthecellsthatneed
aloneRTcodesimplementingAMRexist:someexamplesare to be refined are clustered together to form disjoint rect-
RADAMESH (Cantalupo & Porciani 2011), which is a Monte angular patches, BSAMR schemes are called Patch-based
Carlo (MC) RT code with a ray-tracing scheme, and FTTE (PBAMR; Dai 2010). Both BSAMR and PBAMR involve
(Razoumov&Cardall2005),whichimplementsaschemeto refining a specific region in the grid rather than a single
perform RT on refined grids in the presence of diffuse and
MNRAS000,1–17(2015)
Enabling Radiative Transfer on AMR grids in CRASH 3
cell, and so tend to be used synonymously in AMR-related • Representation of each level as composition of
literature. cell boxes.
Hereafter, we will limit our discussion to the PBAMR In the PBAMR scheme of CHOMBO each level is composed
scheme and look at the details involved for its usage in our by cells which should be seen as minimum units of space
RT code CRASH. withassignedwidth.Cellsareorganizedinrectangulargrids
To enable CRASH to process AMR-based grids we have called “boxes” each of which occupies a unique location in
adoptedtheopen-sourceAMRlibraryCHOMBO(Adamsetal. thelevelanditisdisjointfromtheothers,i.e.acellatapar-
2011). CHOMBO 1 is actively developed at Lawrence Berkeley ticular refinement level can belong to only one box. CHOMBO
NationalLaboratoryandimplementsaPBAMRschemeon represents these boxes by the Box class. Each box occupies
aC++frameworktosolvesystemsofhyperbolic,parabolic a unique location in the 3D space, and its coordinates are
and elliptic partial differential equations. It has been suc- provided by the smallEnd and bigEnd functions of the Box
cessfully used by many gas dynamic codes such as CHARM class.Arefinementlevelisrepresentedasanarrayofdisjoint
(Miniati & Colella 2007), FLASH (Dubey et al. 2014) and boxes (DisjointBoxLayout) which can be traversed with it-
PLUTO (Mignone et al. 2012), and hence is suitable for our eratorclasses(e.g.,theDataIterator class)toaccessasingle
purposes as well. One could, in principle, also post-process box.
the outputs of codes adopting similar BSAMR frame-work Boxes that lie adjacent to a box at the same refine-
as done, for example, by ENZO (Bryan et al. 2014). This ment level are said to be its neighbors and can be accessed
would,however,requireanintermediatesteptoconvertboth throughtheNeighborIterator class.Subsetofboxescanalso
thegridhierarchyandfileformatsinadataformatsuitable begroupedtogetherinanarraycalled“FortranArrayBox”
for CHOMBO. The development of this interface will be taken (corresponding to the class FArrayBox) to allow an easy
upinthefuture;inthispaperwelimittothedatathathas and fast access to the subset data for retrieval or update
beendirectlyproducedusingtheCHOMBOlibraryorthecodes operations.
adoptingit,forexampleCHARM.Test2,inSection5.2,shows • Interactionsbetweendifferentrefinementlevels.
one example of CHARM output post-processed by CRASH-AMR When a coarsening or refinement is requested at some
2. level L, it operates on certain box(es) with a given refine-
Thelibraryisorganisedintoahierarchyofclasses,each mentratior,alsoimplyingthatthebox(es)canhavemulti-
of which provides a specific functionality for incorporating plechildboxes(i.e.atlevelL+1)andcanlieovermultiple
AMR into a stand-alone code with minimum effort, so that parentboxes(i.e.atlevelL−1).TheboxesatlevelL−1or
the software developers only need to focus on implement- L+1whoseintersectionwithaboxatlevelLisnon-empty
ing the physics. Some necessary but routine tasks, associ- are its parent and child boxes, respectively.
atedforexamplewithgridgeneration,management,refine- When physical quantities representing continuous fields
mentandtime-stepping,areautomaticallymanagedbythe are refined in space (i.e. the result is stored across different
library. Hereafter, we describe some of the functionalities levels)thecontinuityoftheirgradientsmustbeensured.For
implemented by CHOMBO that have been extensively used in thisreasontheAMRschemeof CHOMBOadoptsinterpolation
CRASH-AMR.Foreachofthem,theadoptedC++classisalso and averaging methods at the interface of grids. A typical
indicated. exampleofthesesmoothingoperationsoccurwhenthedata
onthegridsareloadedtosetuptheinitialconditions(ICs)
• PBAMR organisation in a hierarchy of levels.
of a hydro simulation. To initialise the finer grids from the
The PBAMR-AMR scheme implemented by the library
existing coarse grids, the FineInterp class is used, while to
consistsofvariousrefinementlevelsLorganizedinahierar-
updatethecoarsegridswiththedataonthefinergridsthe
chy starting from the base level (L=0), and extending up
CoarseAverage class is adopted.
to a final level with index L = N −1, where N is the to-
• Data storage and grid I/O. CRASH-AMR adopts the
talnumberoflevelsinthehierarchy.Eachlevelhasitsown
HDF5 data format standard3 to store the RT results, so
resolution defined by r ·rN, where r is the resolution of
0 0 thattheycanbeeasilypost-processedandvisualisedbyus-
level L=0 and r is the “refinement ratio”, i.e. the ratio of
ingstate-of-the-artvisualisationsoftwarelikeVisit4 orPar-
the resolution between two contiguous levels. While the en-
aview5.
tireAMRschemeisrepresentedbyaC++classAMR,each
refinement level L in the AMR hierarchy is implemented in
CHOMBO by the class AMRLevel. Each AMRLevel class im-
3 RADIATIVE TRANSFER CODE CRASH
plementspointerstoitsrelevantparentinthehierarchy(i.e.
levelL−1)anditschild(i.e.L+1)toallowaneasytraversal CRASHisa3DMCRTcodethatcanself-consistentlyfollow
of the entire hierarchy of refinement levels. the formation and evolution of ionised regions created by
sources present in a static and inhomogeneous gas environ-
ment;thegasconsistsofH,Heandmetals.Thetemperature
1 ChomboisaSwahiliwordmeaning“tool”or“container”.
2 The adoption of output data produced by codes implement- evolutionofthegasiscalculatedself-consistently.Addition-
ing different AMR schemes (e.g. RAMSES, which uses a cell-based ally, the code can account for an arbitrary number of point
scheme (Khokhlov 1998) could require a non-negligible effort in sources, as well as a UV background.
convertingboththegridrepresentationsandthefileformats.The Our work is based on CRASH version 3 (CRASH3; refer
ray-tracingscheme,describedinSection4,shouldalsosensitively
adapt to the different nesting geometry of the AMR levels. The
adoption of different AMR schemes is then beyond the scope of 3 http://www.hdfgroup.org/
thefollowingpaper. 4 https://wci.llnl.gov/codes/visit/home.html
5 http://www.paraview.org/
MNRAS000,1–17(2015)
4 Hariharan N. et al.
to Graziani et al. 2013 and references therein for more de- and its subclasses) to store the physical variables that have
tails), and the developments presented here contribute to a spatial representation, e.g. n ,x ,x ,x and T.
gas HII HeII HeIII
CRASH-AMR by using CRASH3 as baseline but, for simplicity, IntheraytracingalgorithmimplementedinCRASH-AMR
without the inclusion of metals. In this section, we discuss (seeSec.3formoredetails),theinteractionofradiationwith
the CRASH code briefly, and we provide only the details rel- matter is computed in each crossed cell by solving the ioni-
evant to the AMR implementation. sationandtemperatureequations.Thisimpliesthatwhena
CRASHworksbyassigningtheICsontoastatic,regular photon packet propagates through the domain it can cross
3Dgridwhichspecifiesthegasnumberdensityn ,temper- many refined regions containing multiple PBAMR CHOMBO
gas
ature T, the H, He ionisation fractions (x ,x ,x ), patches. Moving through different AMR layers implies that
HII HeII HeIII
the source coordinates, luminosity L and spectral energy the instances of the Box classes representing collection of
distribution(SED)S.Theradiationfromeachsourceisdis- cellsateachrefinementlevelarecontinuouslyaccesseddur-
cretised into photon-packets represented by N frequency ing the travel of each photon packet (see Sec. 3) to update
ν
bins, each containing N photons as determined by the and store the physical quantities.
p,ν
SED,whicharepropagatedalongtherayscastedinrandom The computational cost of a continuous and inefficient
directions from the point sources. The simulation proceeds access to the CHOMBO library could impact the global RT
byemittingphoton-packetsfromallthesourcesandpropa- performance.Infact,dependingonthechosenresolutionin
gating them along the rays until the end of the simulation space and the maximum refinement level provided by the
time. gas dynamics simulation, a single ray could traverse a large
WediscussbrieflytheraytracingroutineofCRASHinthe number of cells spanning different refinement levels. This
following paragraphs. Consider a ray along which a packet makes the box iteration computationally inefficient when
propagates by crossing a series of cells. For each cell l that repeated for the large number of rays required by the MC
is crossed, CRASH calculates the casted path δ and the cor- convergence (typically larger than 107). Note that this is
l
responding optical depth of the cell as not the way PBAMR libraries are normally used in hydro
codestoaccesstheinformation:thepatch-basedschemeim-
plemented in CHOMBO is in fact very efficient in the manage-
τ =τ +τ +τ
HI HeI HeII
(1) ment of memory and parallel computational resources, but
=[σ (ν)n +σ (ν)n +σ (ν)n ]δ,
HI HI HeI HeI HeII HeII l provides information at the grid level instead of at the cell-
where n and σ are the number density and cross section based quantities, as required by the CRASH-AMR RT scheme.
A A
of the absorber A = HI, HeI, HeII. If the packet reaches As further complication, a realistic RT simulation gener-
the cell with a photon content of N , then the number of allyinvolvesanirregulardistributioninspaceofthesources
γ
photons absorbed in the cell is given by fromwhichalargenumberofraysisemittedinrandomdi-
rections, implying that the boxes at each refinement level
are not accessed contiguously. Consequently, the standard
Nl =N (1−e−τ). (2)
γ γ interface provided by the CHOMBO library cannot be simply
Nl is then used to calculate the ionisation, recombination re-used in CRASH-AMR. To resolve this issue we have devel-
γ
fractionsandtemperatureequationsthatregulatethephys- opedanovelFortrandatastructureinCRASH-AMR,minimis-
icalstateofthegas(Section2.4ofMasellietal.2003).The ing the run-time overhead to access and iterate the AMR
angular direction of the ray and coordinates of the current layers, as described below.
cell are used to calculate the coordinates of the next cell During each ray traversal the photon-packet informa-
that the ray will cross; this is repeated until the photon tionhastobepropagatedthroughtheAMRgridhierarchy.
content in the packet is extinguished or, if periodic bound- If we assume that the source lies in the highest refinement
aryconditionsarenotapplied,thepacketexitsthegrid.We level, then the ray will propagate, starting from that level,
finallynotethatsinceweuseCRASH-AMRinapost-processing to the coarser levels below and move back to finer levels if
mode, operations of data smoothing, described in the sec- present. At each step of the RT simulation, the data struc-
tion above, are confined to the initialisation of the RT and ture needs to know the refinement level a cell belongs to,
do not impact the ray tracing algorithm. whether the cell is refined or not, and which box contains
therefinedcell;thewayallthesequantitiesareaccessedand
themappingbetweenCRASHandCHOMBOdataisdescribedin
thefollowingparagraphs.Hereafter,wewillrefertothegrid
4 COUPLING CRASH AND CHOMBO
with the lowest resolution as ‘base grid’, while the refined
In this section we describe how CRASH and CHOMBO are cou- grids will be referred to as ‘refined levels’.
pled to implement CRASH-AMR; we discuss both the adopted Figure 1 shows the data representation in CHOMBO on
methodologyandthesolutionswefoundtothevarioustech- theleft-handside(seealsoSection 2formoredetails),and
nical issues that occurred during the code development. the CRASH equivalent on the right-hand. We also use simi-
CRASH-AMR is implemented in Fortran 2003, while lar colours in both sides to represent corresponding boxes
CHOMBO is a C++ code. We have created a C interface be- at a given refinement level. A simplified but representative
tween Fortran and C++ to allow CHOMBO to communicate CHOMBOhierarchy,consistingofthebaselevel0anditsrefine-
and share information with CRASH-AMR by using the inter- ments from 1 to L−1, is shown in the picture; for clarity
operability features implemented in the programming lan- purposes we just represent the boxes at levels 0, 1 and 2.
guages specifications. As discussed in Section 2, CRASH-AMR To help the reader in connecting this picture with the ab-
uses the grid representation of CHOMBO (i.e. the Box class stract data representation provided in Section 3, we point
out that each refinement level (dashed boxes on the left)
MNRAS000,1–17(2015)
Enabling Radiative Transfer on AMR grids in CRASH 5
CHOMBO This can then be used to index into the arrays in CRASH, to
gettherightboxandthephysicaldatastoredatthesource
coordinates. Then, during the propagation of the photon
Level 0 packet, at each cell crossing, the following scenarios apply:
(base grid)
(a) theraymightescapethegridandthenwenolongerfol-
low it unless periodic boundary conditions are applied;
(b) the photon content of the packet is completely ab-
sorbed and then the propagation stops;
(c) the ray crosses the cell and enters a new cell at the
same AMR level;
(d) theraycrossesthecellandentersanothercellatafiner
(or coarser) AMR level.
While cases (a) and (b) do not need further comments, for
case (c) the new cell might lie in the same box or it might
Figure 1. Interface between CRASH and CHOMBO, the grid hierar- enter a new one. In the former case, we just continue the
chyinCHOMBObeingreflectedthroughthedatastructurebuiltin ray propagation as described earlier; in the latter case, we
CRASH.TheAMRgrids,storedasanarrayofboxesinCHOMBO,are
use the coordinates of the new cell the ray is in and search
storedinCRASHonalevelbasis,withpointerstothegriddataat
the neighbors of the box that the ray was previously in, at
eachlevel.EachboxhasanassociatedlocalID andglobalID.See
thesamerefinementlevel,fortheboxthatcontainsthenew
textformoredetails.
cell.Finally,case(d)needsadifferentapproachbecausethe
cellopticaldepth,calculatedusingthecastedpath,depends
is implemented in computer memory by an instance of the on the refinement level the ray is crossing through (see Eq.
AMRLevel class, while the array of boxes at each level is 1). Here, again, different scenarios apply:
implementedbytheDisjointBoxLayout class.Eachbox,for
example the B(0,0) at base level 0 (red box on the left), is (1) the ray enters a finer level. Given the new cell coordi-
an instance of the Box class. nates, we check if the globalID stored in the cell is the
TheCRASHcounterpartoftheAMRhierarchyismapped same as the globalID of the box the ray was in. If not,
on the right-hand side of Figure 1: the base grid is repre- this indicates that the cell is covered by a refined cell.
sented as refinement level 0, while the L AMR grid levels The globalID stored in the new cell is then used to find
are mapped with an array containing pointers to specific the refined box and the corresponding cell. Finally the
propertiesofeachbox.Theboxesateachlevelareuniquely ray is moved to the finer level;
identified,onbothsides,byanassociatedlocalID andglob- (2) therayentersacoarserlevel.Wesearchtheparentlist,
alID. We first use the localID to get direct access to the of the box the ray was previously in, for the new box
right box in the DisjointBoxLayout array, and then to the containing the new cell.
physicalvariablesofinterestduringtheRTsimulation.The
Once the new box is found in the neighbor list for case
correspondingglobalID alongwiththerefinementlevelpro-
(c)orintheparentlistfor(2)werecursivelymovetheray,as
videsanindexintooneofthearraysontheCRASHside.For
incase(1),toafinerleveliftheneighbororparentisrefined.
example,boxB(0,1)atrefinementLevel1hasalocalID of0
The same procedure is repeated until cases (a) and/or (b)
and a globalID of 1. The value of the localID indicates that
apply, and the photon packet propagation stops.
B(0,1) is at position 0 of the disjoint box array storing all
Figure 2 shows the logical workflow described above.
theboxesatthatlevel,andthephysicaldatacanbefinally
In this diagram the solid lines of the flowchart indicate the
accessedthroughitsFArrayBox (seeSec.2formoredetails
standardstepsoftheCRASHraytracingalgorithm,whilethe
on the classes). The start and end coordinates of the box,
dashed lines indicate the steps performed by the CRASH-to-
which determine its size, are also stored to allow the ray
CHOMBO interface to move across the various CHOMBO boxes
tracing algorithm to recognise if a ray has exited a box at
(CHBox) of the AMR side. The many scenarios described
a given refinement level. Along with the physical data, we
above can be visually followed by colours, as shown in the
store in each cell the globalID of the box it belongs to. For
legend of the figure.
cells that are covered by the cells of a refined box, we store
Itisimportanttonotethat,althoughweuseCHOMBOto
theglobalID oftherefinedbox.Thisisusedtodetermineif
initialise and store the AMR grid, once this data has been
a cell that the ray is passing through is refined or not. Ad-
mapped onto the CRASH side, our implementation does not
ditionally, since a PBAMR scheme allows a refined box to
callanyCHOMBOroutinesduringtheRTsimulation.Thedata
lie over multiple coarse boxes, we keep a list of all parents.
structure is used purely to cross the levels and have a fast
This is done by storing, for each box, the globalIDs of its
access to the grid data. As a consequence of this architec-
parent(s). Finally, we also store a neighbor list containing
tural choice, there is no overhead of using CHOMBO during
the globalIDs of all boxes that are neighbor(s) to a box.
theray-tracingroutinebutthetimeneededtofindtheright
The following paragraphs describe how the new data
cell.
structureisusedduringtheraytracingalgorithm.Firstnote
Asafinalconsideration,wewanttoemphasisethatthe
that the sources emitting photons do not move across the
above features are included as a separate functionality, al-
grid during a RT simulation. Hence, given the refinement
lowingtheusertoenableordisabletheuseofAMRgridsto
level the source lies in and its coordinates we need to look
do RT simulations and use the traditional data storage set-
fortheglobalID ofthebox,containingthesource,onlyonce.
tingupthesimulationICsonlyonthebasegrid.IftheAMR
MNRAS000,1–17(2015)
6 Hariharan N. et al.
Start packet loop
NO
Packetindex <= Np End packet loop
YES
Emit photon-packet
Solve RT equations for current cell
Compute next cell crossing
Has the ray exited the NO
CHBox?
YES
YES
Next cell in neighboring CHBox?
NO
Next cell covered by finer cell?
Next cell in coarser level?
NO
Ray has exited the cosmological
grid, no longer followed
Figure 2. CRASH-AMR flowchart describing the interplay between ray tracing and the CRASH-to-CHOMBO interface when a photon packet
travels through many AMR levels. The boxes with solid lines indicate the CRASH-side of the algorithm, while those with dashed lines
indicatetheCHOMBOinterface-side.Thesolid/dashedlinesconnectingtheboxesindicateiftheCRASH/CHOMBOinterface-sideofthealgorithm,
respectively, is being called. Various colours refer to the different algorithmic scenarios that could happen during a photon packet
propagation.Seetextformoredetails.
functionalityisenabled,thentheusercaneitherrunsimple In Test 2, we apply CRASH-AMR to a realistic density field
tests with pre-defined refinement criteria, or more realistic from the CHARM simulations described in Miniati & Colella
cases with AMR grids provided by hydro codes, as shown (2007). Additionally, in Appendix B we take a further look
inthefollowingsection.Byadoptingpre-definedrefinement at the performance of the code in terms of run times and
criteria, the user could decide, for example, to refine an ar- correctnessofresults.Henceforth,weusedtorepresentthe
bitrary part of the base grid and set up specific test cases, comoving distance from a point source, the units in kpc or
whileforrealisticgasconfigurations,therefinementisgener- Mpc are indicated accordingly.
allydeterminedbythehydrocode,andCRASH-AMRoperates
in post-processing mode.
5.1 Test 1: Str¨omgren sphere in a H+He medium
We have set up a test equivalent to Test 2 of the Radia-
5 TESTS AND RESULTS tiveTransferCodeComparisonProject(RTCCP;Ilievetal.
2006b).Thetestsimulatestheevolutionofanionisedregion
In this section we show the results of some of the tests we
aroundasinglepointsourcelocatedatthegridorigin(1,1,1)
haveperformedtoguaranteethereliabilityofthenewAMR
in a box of side length L = 6.6 kpc and mapped on a grid
implementation. We run a number of test cases in idealised
of 1283 cells. The source is assumed to be steady with an
configurations; in Test 1 we compare CRASH-AMR with the ionising rate of N˙ = 5·1048photons·s−1 and an associ-
AMRfunctionalitydisabledtoCRASH3,andthenwecompare γ
ated black-body spectrum at temperature T = 105 K.
resultswith/withouttheAMRfunctionalityenabled.These BB
The volume is filled by a uniform and static gas of number
set-ups are useful to check the numerical noise introduced
density n = 10−3cm−3, containing H (92% by number)
by the presence of the AMR grids on the RT algorithm. gas
and He (8%). The gas is assumed to be fully neutral and
MNRAS000,1–17(2015)
Enabling Radiative Transfer on AMR grids in CRASH 7
1 1
0.1 0.1
x , no AMR
xHII 0.01 xxHI,I, CCRRAASSHH33 xHII 0.01 xxHHII,I, nAoM ARM (R1 r.l.)
HI HII
0.001 xHII, CRASH-AMR (no AMR) 0.001 xHI, AMR (1 r.l.)
x , CRASH-AMR (no AMR)
HI
0.0001 0.0001
∆[%] 0.01 xxHHIII ∆[%] 101 xxHHIII
0.1
0.001 0.01
4 4
T, CRASH3 T, no AMR
4K] / 10 2 T, CRASH-AMR (no AMR) 4K] / 10 2 T, AMR (1 r.l.)
T [ T [
1 1
1 T 1 T
∆[%] 0.1 ∆[%] 0.1
0.01 0.01
0.001 0.001
1 1
0.1 0.1
x 0.01 xHeII, CRASH3 x 0.01 xHeII, no AMR
xHeIII, CRASH3 xHeIII, no AMR
0.001 xHeII, CRASH-AMR (no AMR) 0.001 xHeII, AMR (1 r.l.)
xHeIII, CRASH-AMR (no AMR) xHeIII, AMR (1 r.l.)
0.0001 0.0001
%] 0.1 xxHHeeIIIII %] 1 xxHHeeIIIII
[ [
∆ 0.01 ∆ 0.1
0.001 0.01
0 1 2 3 4 5 6 0 1 2 3 4 5 6 7
d [kpc] d [kpc]
Figure 3. Spherically-averaged profiles at time t = 500 Myr Figure 4.Spherically-averagedprofilesattimet=500Myrfor
for Test 1a. The colours refer to CRASH3 (green) and CRASH-AMR Test1b.ThecoloursrefertoCRASH-AMR(AMRdisabled;red)and
(AMR disabled; red). The bottom sub-panels show ∆ between CRASH-AMR with one refinement level (1 r.l.; blue). The bottom
theCRASH3andCRASH-AMR(AMRdisabled)results.Top:Profiles sub-panels show ∆ between the CRASH-AMR (AMR disabled) and
ofxHII (dashedlines)andxHI (dottedlines).Middle:Profileof CRASH-AMR(1r.l.)results.Top:ProfilesofxHII(dashedlines)and
T (dash-dot lines). Bottom: Profiles of xHeII (solid lines) and xHI (dashed lines). Middle: Profile of T (dash-dot lines). Bot-
xHeIII (dot-dashlines). tom:ProfilesofxHeII (solidlines)andxHeIII (dot-dashlines).
at a temperature T=100 K, which is then calculated self- showsspherically-averagedphysicalquantitiesasafunction
consistently with the progress of ionisation for a simulation ofd,togetherwiththepercentagedifference(∆)betweenthe
time t = 500 Myr, starting at redshift z = 0.1. We out- CRASH3 and the CRASH-AMR results. We define ∆ = (R −
sim ref
put the results at intermediate times t = 10, 50, 100 and R)·100/R ,whereR andR refertoresultsof CRASH3
i ref ref i
200 Myr as in the original set-up. It should be noted that andCRASH-AMR,respectively.FromtheFigureitisclearthat
simpler tests (e.g., with a gas composed by H only or with thereisnosignificantdifferencebetweenthetwocodes,while
the temperature kept constant) have been run as well, and the spikes that we see in ∆ (with a maximum of ∼1%, but
give results similar to those discussed in the following. mostlybelow0.1%)areduetonumericalartefactscausedby
optimisation of CRASH-AMR involving rearranging of double-
precision floating point arithmetic expressions, and are not
5.1.1 Test 1a: AMR disabled
duetothechangesassociatedwithCHOMBO.Thisshowsthat
ToverifythatthechangesdonetoenableRTonAMRgrids theAMRfeatureinCRASH-AMRisisolatedfromtherestofthe
do not introduce any numerical noise, we have run Test 1 codeandcanbedisabledwithoutintroducinganynumerical
withAMRdisabledinCRASH-AMRandcomparedtheresults noise into the results.
to those from CRASH3.
The outcome is given in Figure 3, where each panel
MNRAS000,1–17(2015)
8 Hariharan N. et al.
5.1.2 Test 1b: AMR enabled of this simulation can be immediately used as an input to
CRASH-AMRbyextractingthenecessaryinformationfromthe
Our next step has been to run the previous test with AMR
HDF5 metadata. As the simulation does not provide infor-
enabled,refining1003 cellsaroundthesourcebyafactorof
mation on the star formation, we define the point source
2,withtheICssetuponthebasegridaswellasonthere-
locations associating them to the gas density peaks at the
finedgrid.Thesourceisnowsurroundedbyacartesiangrid,
most refined level.
at an equivalent resolution of 2563, of ∼5 kpc. The results
We set up the following RT simulations:
are given in Figure 4, where each panel shows spherically-
averagedphysicalquantitiesasafunctionofd,togetherwith (a) multiplepointsourcesplacedatthehighestrefinement
the ∆ between the results from CRASH-AMR (AMR disabled, level, at locations far enough so that when moved to
R ) and CRASH-AMR with one refinement level (1 r.l., R). lower refinement levels they do not gather. The sources
ref i
Here we find that ∆ can be as high as 10% for x and are monochromatic, the gas temperature is kept con-
HI
x cells very close to the source (d∼ 0 - 1 kpc), and at stant throughout the simulation;
HeIII
adistanceof3-6kpc,inthepartiallyionisedregion.More (b) as(a),butnowthepointsourcesareplacedatlocations
typicalvalues,though,donotexceed1%forallthephysical close enough so that they can be gathered at the lower
quantities. refinement levels;
FromSection5.1.1weknowthat,withtheexceptionof (c) same point source locations as (a), but with a black-
thespikesobservedforx ,thedifferencesseenarenotdue body spectrum, and the gas temperature is calculated
HI
to the implementation in CRASH-AMR, but they must rather self-consistently with the progress of ionisation during
be associated to the higher grid resolution in the refined the simulation.
region. Also, these differences are not due to any spurious
In all cases the point sources are located within high
noiseproducedbythisparticularmethodofcalculatingthe
density peaks, chosen to ensure that the criteria mentioned
average, but is rather due to the higher resolution being
abovearesatisfiedforourtests.Sincethebasegridisrefined
used. When we have a grid at a higher refinement level all
onlyinthecentralregionwithinwhichthehigherrefinement
therefinedcellsthatcoveracoarsecellmightnotliewithin
levelsalsolie,theresolutionatlargedistancesfromthepoint
a given radius from the source and thus not contribute to
sources is the same for all cases.
the spherical average. This effect is more prominent in the
Wesetareferenceionisationrate,N˙ ,forthesource
partiallyionisedregionswheretwoadjacentcellsmightnot γ,ref
in the highest gas density peak. For the other sources i:
have the same e.g. x values, unlike a fully ionised region
HII
wherexHII is1.Ifacelldoesnotliewithinthegivenradius, N˙ = N˙γ,ref ·mi, (3)
the average value will differ. Also note that Helium has a γ,i m
ref
recombination rate five times higher than that of H. As a
where m and m are the mass in the cell containing the
result,theheliumcomponentsshowmoresensitivityandwe ref i
reference source and source i, respectively. The initial tem-
get differences of ∼ 1 - 10% in many cells. CRASH-AMR will
perature is T = 100 K and the gas is assumed to be fully
then provide a better description of the regions with fully
neutral. The simulation time is 500 Myr.
ionisedhelium,generallymoreconfinedtothebrightestand
To emphasise the advantage of an AMR scheme, we
x-rays luminous sources (e.g. quasars).
compareresultsofsimulationsrunwithdifferentrefinement
Additional tests showing the dependence of the results
levels.Additionally,asmentionedabove,thesourcesarelo-
on the grid resolution are detailed in Appendix A.
cated at the highest refinement level, so if one or more of
them lie within the same cell at the coarser levels, we con-
5.2 Test 2: a realistic density field siderthemtobeasinglesourcewithluminositygivenbythe
sum of the corresponding luminosities at the finest level.
In this section we apply CRASH-AMR on a density field snap- ToensureagoodconvergenceoftheMCcode,wesam-
shot obtained from a simulation run within the Santa Bar-
pletheradiationfieldwithanumberofphotonpacketshigh
bara Cluster Comparison Project, where the formation of
enoughtoreachconvergenceforeachtestcaserunatdiffer-
a galaxy cluster in a standard CDM universe is followed
entrefinementlevels.WefindthattheMCschemeconverges
(Frenk et al. 1999). The simulation has been performed by with108 photonpacketspersource(0.07%differenceinvol-
the hydro code CHARM (Miniati & Colella 2007) in a box ume averaged x values between two test cases with 108
HII
size L = 64 Mpc (comoving) at redshift z = 0.1. The and 109 photon packets per source). However, the conver-
cosmological parameters are Ω = 1,Ω = 0.1, Ω = 0
m b l gence of the MC scheme is very much problem dependent,
and H =50 km s−1 Mpc−1. The simulation is initialised at
0 hence we do not discuss this further.
z = 40 with a base grid of 643 cells representing a box of
643 Mpc3 comoving, and a grid of 1283 cells placed at the
centre of the base grid and representing a box of comov- 5.2.1 Test 2a: multiple sources set far apart - constant T
ing length 32 Mpc. Only the central region is refined based
Here we place twenty point sources far enough from each
on a local density criterion, with a refinement ratio of 2.
other to remain separate at all refinement levels. This con-
At the end of the simulation there are six refinement lev-
figuration tests the effect of grid resolution on the RT sim-
els in total, along with the base grid. The cell width at the ulation.WeadoptN˙ =8·1053photons·s−1,eachpoint
coarsest level is 1 Mpc, while that at the finest level, with γ,ref
source is monochromatic with E = 13.6 eV, and the tem-
an equivalent resolution of 40963 cells, is 15 kpc. The code ν
perature is kept constant throughout the simulation. For
CHARM adopts the CHOMBO library to implement the AMR
simplicity, we consider a H only gas.
functionality, and the HDF5 files available from the output
Figure 5 shows the maps of x created in simulations
HII
MNRAS000,1–17(2015)
Enabling Radiative Transfer on AMR grids in CRASH 9
Figure 5.MapscutthroughthesimulationvolumeforTest2a.Top:MapsofxHII attimet=100Myr.Second from top:Mapsof
xHII attimet=200Myr.Third from top:MapsofxHII attimet=500Myr.Bottom:Mapsofngas,thedottedlinesrepresentthe
extentofthedifferentrefinementlevelsassociatedwithngas (thebasegridisnotseenhere):magenta(1str.l.),orange(2ndr.l.),cyan
(3rd r.l.), green (4th r.l.), yellow (5th r.l.) and red (6th r.l.). From left to right, the columns refer to simulations run with three, four,
fiveandsixrefinementlevels(seetextformoredetails).
with different refinement levels at times t = 100, 200 and els and differences in the ionisation pattern, as the sources
500Myr.Inthebottompanelswealsoshowthegasnumber are able to maintain their surrounding regions fully ionised
densityfield(n ).Dottedlinesrepresenttheextentofthe against the progressively steeper changes in density intro-
gas
differentrefinementlevels(seethecaptionformoredetails). duced by AMR. On the other hand, the extent of the fully
Forreference,wealsoshowthelocationofthemostluminous and partially ionised HII regions shows obvious differences,
point source. Note that, because the RT is done in post- astheygetsmallerandsharperwithhigherresolution.This
processing,thegasconfigurationdoesnotchangeduringthe is due to the larger changes in density and gas recombina-
ionisation evolution. tion rate (which increases by a factor of 3.5 between 6 and
At t=100 Myr the x maps are very similar. Differ- 3 r.l.) present in the more refined grids. As a result, the
HII
ences become more visible at t = 200 Myr, where separate escape of ionising photons becomes more difficult, delaying
bubblescanbeseenontherightsideoftheboxaswegoto the propagation of the ionisation fronts.
higher refinement levels. The largest differences are present Finally,notethatthepresenceofmultiplepointsources
at the final time t=500 Myr. From a comparison between on different planes of the cube and resolved by differ-
theionisationandgasnumberdensitymaps,wecanobserve ent AMR layers, creates an intricate combination of three-
nodirectcorrelationbetweenpositionsoftherefinementlev- dimensionalRTeffectsinthefinalconfigurationoftheover-
MNRAS000,1–17(2015)
10 Hariharan N. et al.
Figure 6.Mapscutthrough thesimulation volume for Test2a. Top:Mapsof xHI attime t=100Myr.Second from top:Mapsof
xHI at time t=200 Myr. Third from top: Maps of xHI at time t=500 Myr. Bottom: Maps of ngas, the dotted lines represent the
extentofthedifferentrefinementlevelsassociatedwithngas (thebasegridisnotseenhere):magenta(1str.l.),orange(2ndr.l.),cyan
(3rd r.l.), green (4th r.l.), yellow (5th r.l.) and red (6th r.l.). From left to right, the columns refer to simulations run with three, four,
fiveandsixrefinementlevels(seetextformoredetails).
lapping HII fronts. This is more evident in the case with sitioncorrespondingtoalmostneutralgas,whileorange-to-
5 and 6 refinements, in two distinct regions. The HII re- redareasmarkalmostionisedgas(i.e.xHI (cid:46)10−5).Ascom-
gion on the left is in fact formed by a point source lying on mentedabove,theAMRrefinementsclosetothesourcescan
a plane different from the one hosting the most luminuous resolve more over-dense structures. While the radiation is
pointsource,andevolvesdifferentlywithincreasingnumber sufficient to substantially ionise the entire area (x (cid:38)0.9)
HII
ofrefinementlevels.Thisprovidesafinalbubbledistribution and to allow the escape of ionizing photons in far under-
and overlap in space which is very sensitive to the number dense voids (see for example the one in the lower right side
of adopted AMR refinements. of the panels), many inner regions still show an intricate
The effects of resolving more and more gas clumps by pattern of residual neutral gas: a large yellow area preserv-
progressively increasing the AMR resolution are better ap- ing a residual fraction x ∼ 4×10−4 surrounds the red
HI
preciated by showing complementary x =1−x maps. spots and blends into green and cyan areas when the neu-
HI HII
ThisisdoneinFigure6,withthesamepanelorganizationof tral fraction progressively increases up to x ∼ 10−2 and
HI
thepreviousfigure.Firstnotethattheregionshownhereis x ∼ 10−1, respectively. At megapaserc scales the struc-
HI
closertotheimagecenter;thisisdonetobetterzoom-inthe turesinthevariouspanelsdifferforonlyfew,minordetails,
spatialdistributionoftheneutralfractionatalltimes6.Also while the external contours show a clear reduction of the
notethatinthesepanelsthecolourpaletteindicatesneutral ionised gas (e.g. focus on the cyan area connecting the two
gas fraction in logarithmic scale with a cyan-to-white tran- central HII regions) from left (3 r.l.) to right (6 r.l.).
Also the volume averaged ionisation fraction depends
on the refinement levels used, with x =4.94, 4.74, 4.51
6 Aseffectofthere-centering,thedistancescaleinFigure6does and4.39·10−2 forthe3rd,4th,5thandHI6Ithrefinementlev-
notcorrespondtotheoneinthepreviousfigure.
MNRAS000,1–17(2015)