RAY-RAMSES: a code for ray tracing on the fly in N-body simulations

3.7 MB·

by  Alexandre Barreira

RAY-RAMSES:acodeforraytracingontheflyinN-bodysimulations
Alexandre Barreira,1,2,3,∗ Claudio Llinares,1 Sownak Bose,1 and Baojiu Li1
1InstituteforComputationalCosmology,DepartmentofPhysics,DurhamUniversity,DurhamDH13LE,U.K.
2InstituteforParticlePhysicsPhenomenology,DepartmentofPhysics,DurhamUniversity,DurhamDH13LE,U.K.
3Max-Planck-InstitutfürAstrophysik, Karl-Schwarzschild-Str. 1, 85748Garching, Germany
WepresentaraytracingcodetocomputeintegratedcosmologicalobservablesontheflyinAMRN-body
simulations. Unlikeconventionalraytracingtechniques,ourcodetakesfulladvantageofthetimeandspatial
resolutionattainedbytheN-bodysimulationbycomputingtheintegralsalongthelineofsightonacell-by-cell
basis through the AMR simulation grid. Moroever, since it runs on the fly in the N-body run, our code can
produce maps of the desired observables without storing large (or any) amounts of data for post-processing.
We implemented our routines in the RAMSES N-body code and tested the implementation using an example
6
ofweaklensingsimulation. Weanalysebasicstatisticsoflensingconvergencemapsandfindgoodagreement
1
withsemi-analyticalmethods.Theraytracingmethodologypresentedherecanbeusedinseveralcosmological
0
analysissuchasSunyaev-Zel’dovichandintegratedSachs-Wolfeeffectstudiesaswellasmodifiedgravity.Our
2
codecanalsobeusedincross-checksofthemoreconventionalmethods, whichcanbeimportantintestsof
g theorysystematicsinpreparationforupcominglargescalestructuresurveys.
u
A
I. INTRODUCTION recipes to model the effects of baryonic physics and can be
0
used in the generation of mock data sets to calibrate obser-
2
Observations of large scale structure in the Universe have vationalpipelines. UsingN-bodysimulationstopredictinte-
] been playing a crucial role in getting ever tighter con- grated effects along lines of sight that cover a redshift range
O straints on competing theoretical cosmological models. Per- is, ingeneral, notasstraightforwardasgettingpredictionsat
C haps the most classical type of such observations consists fixedredshiftvalues. Forthelatter,itoftensufficestoanalyse
. in mapping the three-dimensional distribution of galaxies thesimulationoutputatagivensnapshotoftheparticledistri-
h
p (which trace with some bias the total matter distribution) bution,whereasfortheformeritisrequiredthatthesimulation
- withspectroscopicsurveys[1–3]. Thesesurveysmeasurethe results are analysed in a continuous range of redshift values.
o baryon acoustic oscillations (BAO) signal [4, 5] and cluster- Asaresultofthat,N-bodymethodsforintegratedobservables
tr inganisotropiesduetoredshiftspacedistortions(RSD)[6–8], areoftensubjecttoanumberofapproximationsthataremore
s which allow to put constraints on the rate at which the Uni- or less valid depending on the exact observable studied. For
a
[ verse expands and the rate at which structure grows in it. A instance, one of the most popular methods for cosmological
complementary approach to traditional galaxy surveys (and weak lensing simulations consists of projecting the continu-
2 which is the focus of this paper) is to infer the large scale ousmatterdistributionalongthelineofsightintoaseriesof
v
structureofmatterbymeasuringitsintegratedeffectonlight lens planes [21–34]. This naturally erases the details of the
2
1 that travels from background sources towards us. These in- time evolution of the fields along the line of sight. Further-
0 clude shifts in the temperature of cosmic microwave back- more, theseprojectionsassumethatthesuperpositionprinci-
2 ground (CMB) photons caused by inverse Compton scatter- ple holds for the lensing effects of gravity, which is not nec-
0 ing with high energy electrons inside galaxy clusters – the essarily true in theories beyond General Relativity that have
1. so-called Sunyaev-Zel’dovich (SZ) effect [9–12], shifts in nonlinearequations. SimulationsoftheISWeffectalsomake
0 thetemperatureofCMBphotonsastheycrosstime-evolving certainsimplifyingassumptionsinthecalculationofthetime
6 gravitationalpotentials–theso-calledintegratedSachs-Wolfe derivativeofthegravitationalpotentials(seee.g.[35]).
1 (ISW) effect [13–16], and magnification and distortions of Agreatdealofeffortisnormallyputintoassessingtheva-
:
v background light sources as their emitted light bends due to lidity of the approximations made in these numerical meth-
i strongandweakgravitationallensingeffects[17–20]. ods,andingeneral,theyseemtoberobustenough. However,
X
As observational surveys of large scale structure keep at- giventheeverhigherprecisionobservationsthatlieahead,it
r taining higher precision, it is important that our theoretical is desirable that the same observables can be computed with
a
understandingoftherelevantphysicalprocesseskeepsevolv- differentmethods,especiallythosewhicharesubjecttofewer
ing as well. This helps in the interpretation of current data, approximations. This can allow for important checks of any
as well as in the planning of future missions. In comput- residualtheorysystematicsthatcouldstillbepresent. More-
ing these theoretical predictions, theorists typically recourse over,currentN-bodymethodstoprobetheclusteringofmatter
to N-body simulation methods since these are currently the alongthelineofsightsometimesrequiresubstantialamounts
mostaccuratewaytopredicttheclusteringofmatteroninter- of data to be stored before it is post-processed to compute
mediateandsmallscales,wherethedensityfluctuationshave thedesiredsignal. Thisprovidesextramotivationtodevelop
become nonlinear. N-body simulations allow also to include newnumericaltechniquesthatarelighterindatastorage,es-
peciallyinlightoftheupcominggenerationofsurveys,which
willrequirelargesimulationsandmockdatasetsforcalibra-
tionpurposes.
∗Electronicaddress:[email protected] This paper is precisely about a numerical method for in-
2
tegrated cosmological observables that goes beyond existing Togiveaconcreteexample,conventionalweaklensingstudies
techniques in the number of approximations made and data usually employ the so-called multiple lens-plane approxima-
storage concerns. Our method, which is based on the origi- tion,inwhichtheobservablesarecalculatedonlyonaseries
nalideaofRefs.[36,37],isdesignedtotraceraysfromsome of planes, onto which the density field has been projected2.
sourceredshifttotheobserver,ontheflyintheN-bodysim- Although one can always perform convergence tests on the
ulation. Our algorithm is implemented in the adaptive mesh number of planes used (e.g. Ref. [32]), some of the detailed
refinement (AMR) RAMSES N-body code [38] and performs informationonthetimeevolutionalongthelineofsightisin
the integrations along the line of sight on a cell-by-cell ba- generallost.
sis, fully exploiting the spatial and time resolution provided Ourmainmotivationtodevelopthecodepresentedinthis
bytheN-bodycode. Moreover,sinceitrunsontheflyinthe paper was to overcome the two above-mentioned aspects.
N-bodysimulation,itcanproducethedesiredmaps(lensing, Namely, weaimedtoimplementanumericalmethodthat(i)
ISW, etc) without having to output the particle snapshots for computestheintegratedobservablesontheflyinthesimula-
post-processing analysis. The goal of this paper is to intro- tion,therebyavoidingtheneedtostorelargeamountsofdata;
ducethecodeimplementationandillustrateitsapplicationin and(ii)takesfulladvantageofthespatialandtimeresolution
cosmologicalweaklensingsimulations. of the N-body run to compute the integrals along the line of
Theoutlineofthispaperisasfollows.WestartinSec.IIby sight. Our numerical implementation is based on the origi-
describinginmoredetailthereasonsthatmotivatedustode- nalideaofRef.[36]forweaklensingsimulations,whichwas
veloptheraytracingcodepresentedinthispaper. SectionIII later optimized in Ref. [37]. In particular, in the latter work,
explains the main aspects of the implementation of the algo- theauthorsrealizedthatinparticle-mesh(PM)N-bodysimu-
rithmintheRAMSEScode. InSec.IV,weexplaintheformal- lations,theintegratedquantitiescanbecomputedanalytically
ismtoperformweaklensingstudieswithourcodeandtestour onacell-by-cellbasisasthesimulationisrunning. Theseau-
implementationinSec.VforafixedGaussianpotential. Sec- thorsimplementedtheirmethodinregulargridPMcodes. In
tionVIisdevotedtoweaklensingcosmologicalsimulations. this paper, we follow a similar approach, but implement the
We analyse our code results for one- and two-point statistics algorithminthepubliclyavailableRAMSEScode[38],which
ofthelensingconvergence,whereweassesstheimpactofN- canachieveafargreaterresolutionduetoitsAMRnature.
body resolution and different integration methods. We also When designing this code, it was also our goal to make it
investigatethelensingsignalarounddarkmatterhaloesinour general enough so that it could be used as a platform to per-
simulations. Finally,wesummariseandanalysetheprospects formstudiesofothertypesofintegratedobservables,andnot
forfuturedevelopmentsandworkinSec.VII. justlensing. Asaresult,eventhoughinthiscodepresentation
paper we illustrate the code operation for lensing, we stress
thatthealgorithmismoregeneralthanthat. Inshort,thecode
II. MOTIVATINGANEWRAYTRACINGCODE wepresentinthispapercancalculateintegralsoftheform
(cid:90)
Ingeneral,numericalstudiesofintegratedobservablesem- I = K(χ)Q(x,y,z)dχ, (1)
ploy the following three step strategy. First, one runs a N-
body simulation for a given cosmology and stores the parti-
where χ is the comoving distance along some ray trajectory,
cledataataspecifiednumberofredshiftvalues. Second,the
K is an integration kernel and Q is any field that can be de-
outputfromthesimulationsisusedtogenerateamocklight-
terminedinsidethesimulationboxatcoordinatesx,y andz.
cone from some observer to a given source redshift1. Here,
The calculation of different observables corresponds to dif-
oneoftenneedstoemploysomeinterpolationschemetocon-
ferent expressions for K and Q. For instance, for thermal
structacontinuousmatterdistributionfromthesimulationre-
and kinetic SZ studies, Q is related to the density-weighted
sultsthatareavailableonlyatafinitenumberofredshiftval-
temperatureandbulkvelocityofelectronsinclusters,respec-
ues. Finally,raysaretracedacrossthelightconetoprobethe
tively; for ISW studies, Q is given by the time derivative of
distribution of matter along the line of sight. This strategy
thelensinggravitationalpotential; andforlensingstudies, Q
hasbeenemployedmostnotablyinweaklensingstudies(see
wouldbeassociatedwithsecondtransversederivativesofthe
e.g. Refs. [21–34] and references therein), but also in ISW
lensing potential (cf. Sec. IV). For the case of lensing, rays
[35,39,40]andSZ[41–45]relatedwork.
may also get their trajectories bent (although we anticipate
One can identify, however, two less appealing aspects of
herethatthisisnotthecaseinthisfirstversionofthecode).
thisstrategy. Thefirstoneispracticalandrelatestothelarge
Aswecommentedabove, thereisalreadyasubstantialbody
amounts of data that are needed to generate lightcones for
of work available in the literature on these topics. The code
post-processing.Thesecondisrelatedtothelossinresolution
we present here provides a different method to compute the
alongthelineofsightthatfollowsfromanalysingalightcone
that has been constructed from a finite number of snapshots.
2See,however,theapproachofRef.[29],inwhichthelensingquantities
areintegratedusingthethree-dimensionaldistributionofthesimulations
1Therearehoweverwaystocomputethelightconeontheflyinthesimula- (withoutprojectionontoplanes),butwhichisstillonlyavailableatafinite
tion(e.g.Refs.[23,25]). numberofredshifts.
3
same quantities which, amongst other things, can be used in fathercell,togetherwithitssoncells,formaso-called"grid"
importantcross-checksofthetraditionalmethods. or "oct" of the refined level. If one of these eight sons gets
We note in passing that the ray-tracing machinery that we further refined, then it will form a grid of the second refined
installedinRAMSESmayalsoserveasastartingpointtode- level, i.e., itssoncellswillhavecellsizeh/4. Thisseriesof
velop a code that could be applied in radiative transfer stud- gridsaccrossrefinementlevelsisorganizedinatreestructure.
ies(seee.g.Refs.[46–50]andreferencestherein). Thiswill, RAMSES stores the grid and cell IDs in separate arrays and
however, require some modifications to the code presented in a way that (i) given a son cell’s relative position inside its
here, which is primarily oriented for integrated observables grid and the ID of that grid, one can find the ID of the son
alonglinesofsight. cell, and vice versa; and (ii) given the parent cell’s ID one
can find the grid ID. Each level of refinement is labelled by
l. The domain level l = ldomain is defined by 2ldomain =
III. CODEDESCRIPTION N1/3. For instance, if the simulation has N = 5123, then
p p
l =9.Thefirstrefinementlevelislabelledbyl +
domain domain
Inthissection,wedescribethemainpartsoftheraytracing 1 and so on. Another characteristic of the RAMSES code is
code. We start by presenting a quick review of the default that,atrefinementboundaries,thecoarseandfinesidesdiffer
RAMSES code, which is followed by an overview of the ray only by one level of refinement. The size of the time steps
tracingalgorithmandhowitisimplementedinRAMSES.We is determined independently for each refinement level, with
thenexplainwithmoredetaileachofthemainpartsoftheray higher refinements taking smaller time steps. For example,
tracingcode. one of the criteria to determine the size of the time steps is
that the particles should move only by a fraction of the cell
sizetheyarecurrentlyin(cf.Sec.2.4ofRef.[38]).
A. NotationandbasicsoftheRAMSEScode RAMSES is also efficiently parallelised using MPI. When
run in parallel, each grid is "owned" by the same CPU that
Ourraytracingmodulesareinstalledinthepubliclyavail- "owns"theparentcellofthegrid. "Toown"heremeansthat
able AMR RAMSES N-body code [38]. RAMSES employs a theCPUknowsallnecessaryinformationofacell/grid,andis
numberofsimulationparticlesN whichactasdiscretetrac- responsibleforcalculatingalltherelevantquantitiesinsidethe
p
ersoftheunderlyingmatterfield. Thesimulationboxiscov- cell(density,potential,etc.),aswellasthesoncellpositionin-
ered by a three dimensional mesh, on which the density val- side the grid. As domain decomposition strategies, RAMSES
uesarecalculatedusingthecloud-in-cell(CIC)interpolation can employ Peano-Hilbert, angular and planar schemes. Of
scheme given the particle distribution at any time step. The these,thePeano-Hilbertspacefillingcurveisthemostoptimal
code then solves for the gravitational potential field on the forstandardN-bodyruns. However,giventheangulargeom-
mesh,whichcanbefinite-differencedtofindthecorrespond- etryoftheraytracingoperationsthatweaimtoperform,itis
inggravitationalforces. Theforceattheparticlepositionsis beneficial to consider the angular scheme since it distributes
obtainedbyinterpolatingbackfromthemeshusingthesame better the rays accross CPUs (specially when the number of
CICschemetoensuremomentumconservation. Thisisthen raysbecomeslarge).
usedtoupdatetheparticle’svelocityandpositionatthenext Our modifications to the RAMSES code are mostly in the
time step. The whole process is repeated from some initial
form of additional independent numerical modules, which
time(typicallyredshiftz =50−100)toalatertime(usually
keep the base code unchanged (except only for a few inter-
z = 0). Thecubiccellsofthe3Dmeshcangetrefinedifthe faces).WereferthereadertotheRAMSEScodepaper[38]for
effectivenumberofparticlescontainedinthemexceedssome
moredetailsaboutitsoperation.
pre-specifiedthreshold,N . Conversely,thecellsarealso
refine
de-refinedifthenumberofparticlesdropsbelowthatthresh-
old. This AMR nature of the code is useful in cosmological
simulations,becauseitallowsforhighresolutioninregionsof
B. Tiledsimulationboxes
highmatterdensity,whilstsavingcomputationalresourcesin
regionsoflowerdensitywheretheresolutioncanbelower.
The term "domain level" refers to the coarsest mesh that High-resolutionN-bodysimulationsofboxesthatarelarge
regularly covers the whole simulation volume. In RAMSES, enoughtocontainthedistancetravelledbyphotonsfromz (cid:38)
the domain level contains N1/3 cells in each direction3. If a 1 (χ (cid:38) 2 Gpc/h) typically require massive computational
p
cellofthedomainlevelgetsrefined,thenitiscalleda"father resources (see e.g. Refs. [24, 25]). To circumvent this, one
cell"witheightcubic"soncells". Ifthecellsizeofthefather can "tile" together a number of simulation boxes in order to
cell is h, then each of the eight sons has cell size h/2. The fit the whole light bundle [21, 36, 37]. Figure 1 shows an
exampleofapossibletilingscheme. Theobserverliesinthe
boxthatwerefertoasthe"lastbox",asopposedtothe"first
box", which contains the ray sources. The source redshift is
3Inthispaper,wealwaysassumethreedimensionalsystems,althoughsome zs = 1 in this example. Each simulation box takes as input
timesweshallusetwo-dimensionaldiagramstofacilitatetheillustrations the position of the observer w.r.t. its origin, xobs,yobs,zobs.
andexplanations. For example, for the case illustrated in Fig. 1, the observer
4
FIG.1.Exampleofatilingschemeforraytracing.Thex-axispointsintotheplaneofthefigure.TheboxeshavesizeL=512Mpc/h.The
thinredlinesillustratethetrajectoryoftheraysinalightbundlewithopeningangleΥ = 10degfromz = 1. Thedifferentboxesshould
s
alsosimulatedifferentrealizationsoftheinitialdensityfieldtoavoidrayscrossingthesamestructuresatdifferenttimes.Iftheraytrajectories
arestraightlines,thentheboxesinthetilecanbesimulatedsimultaneously.Inthecaseofbendingrays,eachboxcanstartwhentheprevious
box(higherredshift)hasfinishedthecalculations. Althoughitisnotthecaseinthefigure,closertotheobserverwherethebundlecoversa
smallervolume,thesimulationboxescanbemadesmallertoincreasetheparticleresolution.
is located at the center of x-y face (z = 0) of the last box4, boxeslocatedclosertotheobservercanstarttracingtherays
and so we have x = y = 256 Mpc/h, z = 0 for afterreadingtheirpositionsfromhigher-redshiftboxes.
obs obs obs
thatbox. Forthefirstbox,ontheotherhand,thesewouldbe Wenotethateachsimulationboxshouldalsostartfromdif-
x = y = 256 Mpc/h, z = −2048 Mpc/h. Given ferent statistical realizations of the initial density field. This
obs obs obs
the geometry of the light bundle, the ray positions are more way,oneensuresthattheraysdonotseethesamestructures
easilydescribedusingasphericalcoordinatesystemwiththe throughouttheirtrajectories(duetotheperiodicconditionsof
observeratitsorigin thesimulationbox). Finally,althoughnotdepictedinFig.1,
itisalsoworthmentioningthatclosertotheobserver,where
x =χsinθcosϕ,
ray the spatial volume covered by the light bundle is smaller,
yray =χsinθsinϕ, the boxes in the tile can be made smaller to gain resolution
z =χcosθ, (2) withoutsacrificingcomputationalefficiency(althoughnottoo
ray
smalltostillallowlargeenoughstructurestoform).
whereθ ∈ [0,π],ϕ ∈ [0,2π]arethetwoangularcoordinates
on the sky and χ is the radial coordinate. If the rays follow
straight trajectories, then χ(z) is equal to the comoving dis-
(cid:82)z C. Outlineofthecode
tance D (z) = c dz/H(z), with H(z) being the Hubble
c 0
expansionrate,ztheredshiftandcthespeedoflight.
Figure 2 shows a sketch of the flow of calculations in the
Inthetilingscheme, arayisonlytracedinagivenboxin
code.Thefirstoperationoftheraytracingcalculationconsists
theredshiftintervalduringwhichtheraypositionlieswithin
of the initialization of the ray data structure (cf. Sec. IIID).
that box. For example, the integration of the rays in the first
The goal is to identify the ID of the grid that a given ray
boxwouldstartatz =1anditwilllastuntilz ≈0.86,which
belongs to, i.e., determining the physical location of the ray
is approximately when the rays "touch" the face of the box.
within the grid structure (cf. Sec. IIIE). This is performed
Followingthesamereasoning,thesecondboxwouldstartthe
only when the rays start the integration because, as the rays
integrationsatz ≈ 0.86,whichwillcontinueuntilz ≈ 0.60;
movethroughthe mesh, itispossibletodetermine theIDof
andsoonandsoforth,untiltheraysreachtheobserveratz =
thenextcrossedgrid,bysearchingforneighbouringgrids.
0. Naturally,rayslocatedintheoutermostregionsofthelight
Aftertherayshavebeeninitialized,theyaremovedacross
bundlemovefromoneboxtotheotherbeforethemorecentral
the mesh on a cell-by-cell basis (cf. Sec. IIIF), integrating a
rays.Theconditionsforthestartandendofintegrationineach
givenquantityalongthepathinsideeachcell. Asweexplain
oftheboxesareexplainedwithmoredetailinSec.IIIF3. We
inSec.IIIG,theintegrationcanbedoneanalyticallybyusing
notealsothatforstraightraycases,theboxesinthetilecanbe
thevaluesofthedesiredquantityateachcrossedcellcentreor
runsimultaneously,sincetheyall"know"aprioritheposition
atitsvertices. Thelatterhavetobeobtainedbyinterpolation
oftheraysatalltimes. Ontheotherhand,ifraysbend,then
from the cell centres, which is where RAMSES evaluates all
fields(density,potential,etc.) bydefault(cf.Sec.IIIH).
Inagiventimestep,eachCPUmovestheraysthatarecur-
4Weusethesameletter,z,todenoteredshiftandoneofthecartesiancoor- rentlywithinitsspatialdomainuntiloneofthefollowingpos-
dinates.Themeaningofzshouldbetakenbythecontext. sibilitieshappens:
5
FIG.2. Sketchofthecodeflow. TheraytracingroutinesarecalledafterdefaultRAMSEScomputesthefieldvaluesandthenextscalefactor,
butbeforetheparticlesaremoved. Theraytracingroutinesareinitializedonceeachsimulationboxinthetilereachestheredshiftatwhich
theraysstarttheirpropagationthere. Ineachtimestep,eachCPUintegrateseachrayonacell-by-cellbasis(i)untilittravelsthemaximum
alloweddistancelightcantravel,inwhichcasetheCPUgoesdirectlytothenextray;(ii)untilitreachesaCPUdomainboundary,inwhich
casetherayismarkedforcommunicationforanotherCPUtocontinueitsintegration;or(iii)untilitreachestheobserverand/ortheendofthe
box,inwhichcasetherayintegrationintheboxismarkedasfinishedforthatray.
i theraystravelthedistancethatlightcantravelinthattime theraytracingroutinesagainatthenexttimestep. Ourmod-
step; ifications consisted therefore in the development of indepen-
dent modules that do not impact in any way the standard N-
ii theraysreachtheobserver/faceofthebox;
bodypart.
In the remainder of this section, we explain in more de-
iii theraysreachtheendoftheCPU’sspatialdomain.
tail each of the steps and concepts involved in the propaga-
tion and integration of the rays across the mesh. The reader
If (i) is satisfied, then the CPU simply moves on to the next
who wishes to skip these details can jump to Sec. IV, from
ray. If(ii)happens,thentheray’sintegrationismarkedasfin-
whereon we present tests and results from weak lensing ray
ished, and the CPU also proceeds to the next ray. Finally, if
tracingsimulations.
(iii) happens, then the ray is marked for communication and
the CPU still carries on to its next ray. Once each CPU has
dealt with its initial number of rays, it checks whether rays
from other CPUs have been marked to enter its domain, and D. Raydatastructure
whetheritsownrayshavebeenmarkedtoleave. Ifthereare
rays entering and/or leaving the CPU’s domain, then the rel- To implement our ray tracing algorithm in RAMSES we
evant CPUs exchange ray data via MPI communication and need to establish a data structure that links ray and grid IDs
repeattheabovecalculationsfortheincomingrays. Thispro- 5. Onecanthinkofatleasttwowaystodoso. Wecallonethe
cessisrepeateduntilallrayssatisfy(i)or(ii).
Our ray integration routines require as input the field val-
ues given by RAMSES at a given time step, and hence, they
arecalledafterRAMSEScomputesthesequantities. Oncethe 5OncethelinkbetweengridandrayIDsisset,thenthelinkbetweentheray
raytracingcalculationsforthistimesteparefinished,thecode andtherelevantcellismadebycheckinginwhichoctantofthegridthe
proceedswiththestandardN-bodypart,untilitistimetocall rayisin.
6
all the code needs to do is to (i) compute the integral for the
crossed cell and accumulate it in the n-th entry of the quan-
tity’s list; and (ii) update the n-th entry of the grid list with
thegridIDthatcontainsthenewcell(whichcanbethesame
grid). Therayapproachismoreefficientthantheparticleap-
proach, even though it involved developing new strategies to
moveraysandcommunicatebetweenCPUs(whenraysmove
outside a CPU’s spatial domain or when RAMSES does load
balancingtoreassigndomainstoCPUs).
E. Rayinitialization
FIG.3. Sketchoftwopossibledatastructureschemestolinkglobal
rayIDswithlocalgridIDs. Theparticleapproachtreatsraysasa
Thegoaloftherayinitializationistosetupthedatastruc-
differentparticletypeinRAMSESassociatingeachgridwithallthe
rayscontainedit.Thisapproach,whichisbasedonRAMSES’slinked ture described in the previous subsection. Just to guide the
lists, is, however, computationally expensive because ray particles discussion in this section, we assume that the central ray of
travelatthespeedoflight,whichrequiresthelinkedliststobeupdate thebundletravelsinthenegativezdirection.Thelightbundle
toomanytimes. Intherayapproach,twoorderedlistslinkeachray isspecifiedby(i)itsopeningangleinthexandy directions,
tothegriditiscurrentlyin. Ifarayleavesitscurrentgrid,thenall Υ , Υ , respectively; and (ii) the number of rays in each of
x y
thereistodoistoupdatetheentryofthegridlistthatcorrespondsto thesetwodirections,N andN . Theangularpositionsofthe
x y
thatray. raysinthefieldofviewareassignedas
Υ Υ
θ = a (i −1)− a, a=x,y, (3)
a N −1 a 2
"particle approach" and the other the "ray approach". These a
twoapproachesaresketchedinFig.3. andwedefinetheglobalrayIDnumberas6
Incrudeterms,theparticleapproachdetermineswhichrays
areineachgrid(cf.left-handsideofFig.3).Theadvantageof rayid=(iy−1)Nx+ix, (4)
itisthatitenablesdirectuseoftheexistingRAMSESstructure
with i ∈ [1,N ] and i ∈ [1,N ]. Once initialized, the
x x y y
for other types of "particles" (dark matter, stars, sinks, etc.),
globalrayIDstaysthesamethroughoutthesimulation,even
thereforemakingthecodingeasier. Inthisapproach,alinked
whenrayschangeCPUs. Attheendofthesimulation,given
list data structure is used in RAMSES to find the global IDs
arayID,onecanreversetheaboveequationstofindtheray’s
of the particles that lie within a grid, given the grid ID. The
positioninthe2Draytracingmap.
communicationofraysbetweenCPUswouldalsofollowthe
In the RAMSES structure, it is straightforward to retrieve
strategyalreadysetupforothertypesofparticlesinRAMSES.
the spatial location of a grid given its ID, but not the inverse
However,inthecode,onecriterionfordeterminingthesizeof
operation:tofindthegridIDgivenacertain(theray’s)spatial
thetimestepensuresthatdarkmatterparticlesmoveonlybya
location. Toachievethis,onecanloopoverallgrids,andfor
fractionofthecurrentgridsize.Duringthisparticletimestep,
each grid, loop over all rays to identify those that are inside
photons, which travel at the speed of light, can cross many
it. We note that such a "brute force" search may not lead to
grids. Thismeansthatonehastoeitherupdatethelinkedlists
huge overheads to the overall performance of the code since
eachtimerayschangegrids,whichinvolvesalargenumbers
the initialization is performed only once per box. Neverthe-
of operations and memory allocations/deallocations, or dras-
less,wehavedevelopedamoreefficientalgorithmthatloops
tically reduce the particle time step so that rays do not cross
overallgrids,butforeach,onlyloopsoverasmallertargeted
morethanonegridinatimestep,whichwouldmakethecode
rangeofrayIDs.ThedetailsofthealgorithmaregiveninAp-
prohibitivelyslow.
pendixA,butinshort, theideaistocomputethesolidangle
The above drawbacks motivated us to implement the ray thatsubtendsaspherethatcontainsthegrid, withwhichone
approach, in which one determines which grid a ray is in can determine a range of θ and θ . This then enables us to
x y
(cf.right-handsidepartofFig.3). Inthiscase,thedatastruc- significantly narrow down the ray IDs to loop for each grid
ture consists of two ordered lists of global ray IDs and their usingEqs.(3)and(4).
corresponding local grid IDs. Ordered here is in the sense Fornon-firstboxesinthetile,theraysareinitializedatthe
that the global ray ID that is in the n-th entry of the ray list box face, and as a result, if the spatial location of the grid
physically resides inside the grid whose local grid ID is that (whichcanbedeterminedstraightforwardlyinRAMSES)does
ofthen-thentryofthegridlist. Notethat,naturally,ifaray notlieontheboxface,thentheloopoverrayscanbeskipped.
is inside a refined grid, then it also physically lies inside its
parent grid. In the data structure, rays are always linked to
the grid that is at the highest refinement level. With the ray
approach, theintegratedquantitiesareassociatedtoeachray 6Thereisanabuseofnotationwiththesubscriptsxandy. Here,theyde-
IDinthesamewayasgridIDs. Forexample,iftherayinthe notethetwodirectionsperpendiculartothelineofsightandshouldnotbe
n-thentryoftheraylistcrossesonecelltoenteranother,then confusedwiththe3DCartesiancoordinates.
7
theequation7
(cid:0) (cid:1)(cid:0) (cid:1)
(cid:126)r = χ +λ sinθcosϕ,sinθsinϕ,cosθ , (5)
A
whereχ isthecomovingdistanceoftheraytotheobserver
A
atthebeginningofthetrajectoryinthecell(cf.Fig.4).Hence,
(cid:0) (cid:1)
(cid:126)r =χ sinθcosϕ,sinθsinϕ,cosθ (6)
A A
isthepositionvectoroftherayatthebeginningofitstrajec-
tory. Ofallofthesixcellfacestheraycancross,threecanbe
ruled out by the signs of cosθ, sinϕ and cosϕ. For instance,
ifcosθ > 0,thentherayistravellinginthenegativez direc-
tion (θ ∈ [0,π/2]). This means that the ray cannot enter the
neighbouring cell that lies in the positive z direction. Simi-
larly, if sinϕ > 0 (cosϕ > 0), then the ray cannot enter the
neighbouringcellthatliesinthepositivey(x)direction.More
compactly,thefacesfromwhichtheraycanleavethecelllie
ononeofthreeplanes,eachcharacterizedby
x=x =x −sign(cosϕ)h/2,
target cell
FIG.4. Exampleofaraytrajectoryinacell. PointsAandB are, y =y =y −sign(sinϕ)h/2,
target cell
respectively,theendingandstartingpointsofthetrajectory. Inthis
z =z =z −sign(cosθ)h/2, (7)
illustration, these two points are the intersection points of the ray target cell
trajectoryandthecellfaces, butingeneralpointsAandB canlie where h is the cell size and x , y and z are the cell
cell cell cell
anywhereinsidethecell. centercoordinates.Thetrajectoryoftheray,Eq.(5),intersects
eachofthesethreeplanes,respectively,at
(cid:18) (cid:19)
sinϕ cosθ
Note also that since rays are linked always to the finest grid (cid:126)rB,x =xtarget 1,cosϕ,sinθcosϕ
(cf.Sec.IIID),onecanalsoskiptheloopovergridsthatare
(cid:18) (cid:19)
cosϕ cosθ
notatthehighestrefinementlevel. (cid:126)r =y ,1,
B,y target sinϕ sinθsinϕ
As we commented above already, the initialization needs
(cid:18) (cid:19)
onlytobeperformedonceperbox. Duringthecourseofthe sinθcosϕ sinθsinϕ
(cid:126)r =z , ,1 , (8)
ray integrations, when a ray leaves a grid to enter another, B,z target cosθ cosθ
there are other ways to efficiently determine the IDs of the
wherethesubscript denotestheendoftheraytrajectoryin-
B
newgrid.
sidethecell(cf.Fig.4). Thefacefromwhichtherayleaves
In this paper, we limit ourselves to analysing the code re- thecellisthattowhichtherayneedstotraveltheleast.There-
sultsforsmallopeningangles,forwhichtheflat-skyapprox- fore,onecomputes
imationholds. Forfull-skyraytracing,onemaybenefitfrom x
using techniques such as HEALPix [51] to describe the ray Dx =(cid:107)(cid:126)rA−(cid:126)rB,x(cid:107)2 =χA− sintθacrgoestϕ,
distributionacrossthesphere. Weleavethegeneralizationto y
full-skycasesforfuturework. Dy =(cid:107)(cid:126)rA−(cid:126)rB,y(cid:107)2 =χA− sintθarsginetϕ,
z
D =(cid:107)(cid:126)r −(cid:126)r (cid:107)2 =χ − target, (9)
z A B,z A cosθ
and the value of min{D ,D ,D } specifies the target face.
F. Movingrays x y z
Forexample,ifmin{D ,D ,D }=D ,thentherayleaves
x y z x
fromthefacex=x (andanalogouslyforyandz).
target
In this section, we describe how the algorithm determines
Once the exiting face is found, then one can identify the
thepathofarayinacellandhowtheraysaremovedthrough
nextcrossedcellbysearchingforneighboursusingthedefault
themeshacrossparticletimestepsandCPUdomains.
datastructureinRAMSES.Ifthenewcellbelongstothesame
grid, then one needs only to update the cell information. If
ontheotherhandthenextcrossedcellbelongstoadifferent
grid,thenoneupdatesthegridinformationaswell. Notethat
1. Trajectoryinsidecells
thenewgridcanlieatadifferentrefinementleveland/orina
different CPU domain, and this information is also recorded
Giventhedirectionofarayinsideacell,thecalculationof inourraydatastructure.
itstrajectoryandthedeterminationofthefacefromwhichthe
rayleavesthecell(whichdeterminesthenextcrossedcell)is
mostly a geometrical exercise. Denoting by θ and ϕ the two 7Forstraightrays,θandϕcoincidewiththeanglesthatspecifythespherical
spherical coordinate angles that specify the direction of the coordinates of the ray (cf. Eq. (2)). We leave the generalization of our
ray, then the position of the ray can be parametrized by λ in algorithmtofollownon-straighttrajectoriesforfuturework.
8
coming rays. In the case sketched in Fig.5, CPU 2 receives
rayAandintegratesituntilitismarkedforcommunicationto
CPU 3, and CPU 3 receives ray B and integrates it until the
end of the time step. These series of CPU communications
occuruntilallrayshavereachedtheendofthetimestep8. In
ourparticularexample,CPU3receivesrayAandintegratesit
untiltheendofthetimestep.
The way that rays move in the third coarse step is analo-
goustothatoftheothertwo,exceptthatthedistancethatthe
raystravelbeforeRAMSESupdatesthefieldvaluesissmaller,
because of the smaller particle time step. In the current im-
plementationofthecode,allraystravelbythedistancedeter-
minedbytheparticletimesteponthefinestlevel,regardless
of which level they belong to. In principle, this is not nec-
essarysince, ifarayonlycrossesunrefinedcells, thenitcan
bemovedbythe(larger)distanceallowedinthecoarsestep,
if the field values at the corresponding cell are not updated
duringthefinerstepstakenbythecode. Theimplementation
of this in RAMSES is, however, slightly more involved and
therefore we opted to have all rays moving by the distance
determinedbythefinestcellsofthemesh.
FIG. 5. Two-dimensional sketch of two rays moving through the 3. Fromboxtoboxinthetile
mesh. ThecolorsindicatedifferentCPUdomainsandthereddots
denotethestartingandendingpointsofarayinsideacell. Therays Asthebundleapproachesthefaceofnon-lastboxes,some
movefromtoptobottominthefigure. Therearethreecoarsetime
oftherayswillcrossthefaceearlierthanotherrays. Forin-
stepsshown,andthemeshgetsrefinedfromthesecondtothethird.
stance, theoutermostraysinthebundlearethefirsttoreach
Notethatthedistancethatraystravelinatimestepbecomessmaller
the box face, whereas the center rays are the last. To exem-
whentherearerefinements.
plifyhowthecodedealswiththistransition,considerthetra-
jectoryoftheoutermostandcenterraysastheyleaveagiven
2. FromcelltocellacrossdifferenttimestepsandCPUdomains box (call it Box 2) to enter another (call it Box 1), and de-
notebyzoutermostandzcentertheredshiftvaluesatwhichthe
end,2 end,2
outermost and center rays cross the face of the box. In this
Figure 5 shows a two-dimensional sketch of rays moving
case, Box 2 propagates the outermost ray until it reaches its
throughthemesh. Shownarethetrajectoriesoftworaysdur-
face. At this point, the integration for this ray stops, but the
ingthreecoarsetimesteps. Duringthefirsttwothereareno
calculationforthecenterrayisstillongoing. Thesimulation
refinements(left),andinthelastcoarsestepsomecellsgetre-
ofBox2shouldthereforeonlybestoppedforz ≤zcenter. As
fined(right). ThedifferentcolorsindicatedifferentCPUdo- end,2
forBox1,itinitializestheraydatastructureatz =zoutermost
mains. Thereddotsalongtheraytrajectoriesindicatewhere end,2
forallrays,includingthecenteronewhosepositionisknown
therayscrosscellsorreachtheendofatimestep.
beforehand 9. However, it starts by integrating only the out-
In the first coarse time step, ray A starts from the middle
ermost ray. The center ray remains "initialized" at the box
of a cell and is moved to the face of this cell. During this
face, and it only starts being integrated when Box 1 reaches
trajectory, the code integrates the desired quantity using the
z =zcenter <zoutermost.
algorithm that shall be described in Sec. IIIG. Ray A then end,2 end,2
continuesitstrajectoryinthenextcrossedcell, butstopsbe-
forereachingthenewcell’sface,atthepointwheretherayhas
G. Integrationincells
travelledthemaximallyalloweddistanceinthattimestep. In
the meantime, Ray B is moved similarly to ray A, but in the
One of the key parts of our algorithm is the computation
domainofCPU2.
oftheintegralofagivenfieldalongtheraytrajectories. Ex-
Inthesecondcoarsetimestep, CPU1movesrayAtothe
face of the cell and marks it for communication to CPU 2,
which is found to own the next crossed cell (CPU 1 would
thencontinuetomovetheotherraysinitsdomain,ifany). In 8Wenoteinpassingthat,inadditiontothecommunicationsthatoccurwhen
raysleaveCPUdomains,thereareothersituationsthatrequirecommuni-
themeantime,CPU2movesrayBuntilitismarkedforcom-
catingtherayinformation,namelywhenthecodeperformsloadbalancing
municationtoCPU3. AfterCPU2hasdealtwithrayB(and
(redistributionofthespatialdomainsacrossCPUs).
alltheotherraysinitsdomain),itcheckswhetherotherCPUs 9Asmentionedpreviously,thegeometryisfixedforstraightraysimulations
have marked rays to be sent to it and moves these incoming andthereforeBox1"knows"apriorithepositionofallraysatanyredshift.
raysasdescribedabove. Similarly,CPU3alsochecksforin- InthiscaseBox1andBox2canberunsimultaneously.
9
plicitly, we wish to evaluate integrals of the form of Eq. (1), inwhichQ denotethevaluesofQ(x,y,z)atthecellvertices
i
whichwerepeathereforeasyreference, P , as indicated in Fig. 4. In Eq. (14), ∆ , ∆ and ∆ are
i x y z
(cid:90) χs givenby
I = K(χ)Q(x,y,z)dχ, (10)
x −x a+(χ−χ )sinθcosϕ
0 ∆ = ray 1 = A ,
x h h
whereχ =χ(z )(withz thesourceredshift),K(χ)issome
s s s
y −y b+(χ−χ )sinθsinϕ
weightingfunctionorkernelthatcanbeapolynomialofχand ∆ = ray 1 = A ,
Q(x,y,z) is any quantity that can be evaluated in the cells. y h h
ComparedwithEq.(1),Eq.(10)ismorespecificasitencodes ∆ = zray−z1 = c+(χ−χA)cosθ, (16)
theinformationthatwewishtofollowraysfromsomedistant z h h
sourcetoanobserveratredshiftzero. TheintegralofEq.(10)
where h is the cell size, and (x ,y ,z ), (x ,y ,z )
ray ray ray 1 1 1
canbesplitintothecontributionfromeachcellcrossedbya
and(a,b,c)are,respectively,thecoordinatesofthegivenray
ray
(cf.Eq.(2)),thepointP ,andthepointAw.r.t.pointP .
1 1
(cid:88) UsingEqs.(14),(15)and(16),itispossibletocompactify
I = I , (11)
c
theexpressionforQas
cells
inwhich (cid:90) χA Q(χ,θ,ϕ)= (cid:88)4 dN(χ−χA)N−1, (17)
Ic = K(χ)Q(x,y,z)dχ, (12) N=1
χB where the coefficients d depend on a,b,c,α ,θ and ϕ, but
N i
with χA and χB being, respectively, the radial coordinate of notonχ(notealsothatwearenowspecifyingQasafunction
the starting and ending points of the ray trajectory in a cell of spherical coordinates). Their expression is given in Ap-
as the ray travels towards the observer (χA > χB) 10 (see pendix B. Since the dN’s do not depend on χ, one can com-
Fig.4). Toperformtheintegralweneedtobeabletoevaluate bineEqs.(12)and(17)towrite
thequantityQatanygivenpointinsidethecells.Thesimplest
waytodosoistotakethefieldtobeconstantinsidethecell I = (cid:88)4 d (cid:90) χAK(χ)(χ−χ )N−1dχ. (18)
and equal to the value at its center Q , as given by default c N A
c
RAMSES. In this case, Q can be taken out of the integral in N=1 χB
Eq.(12)andthecontributionofeachcelltothetotalintegral Theintegrandintheaboveequationisapolynomialinχ,and
becomes sotheintegralcanbesolvedanalytically[37]. Therefore,the
(cid:90) χA calculation of the integral in the cell reduces to the geomet-
IcNGP =Qc K(χ)dχ, (13) rical exercise of determining the position of points A and B
χB and the evaluation of the quantity Q at the cell vertices. In
wherethesuperscriptNGPstandsfornearestgridpoint. We Secs. IVB and IVD we shall see a concrete example of the
note that, in this way, the fields being integrated do not vary useoftheseformulaewhenweapplyittolensing.
continuouslywhencrossingcellboundaries. Equations(13)and(18)providetwopossiblewaystocom-
Another, more involved, way to construct the field at an putethefinalresult. Thelatterhastheadvantageofallowing
arbitrary point inside the cell is via trilinear interpolation of for the integration of a continuous field when crossing cell
thevaluesofQateachcellvertex: boundaries, by appropriately evaluating the fields at the ver-
tices of the cells (see next section). On the other hand, the
Q(x,y,z)=α +α ∆ +α ∆ +α ∆ +α ∆ ∆
1 2 x 3 y 4 z 5 x y
use of Eq. (13) does not involve evaluating the fields at the
+α ∆ ∆ +α ∆ ∆ +α ∆ ∆ ∆ , (14)
6 y z 7 x z 8 x y z cellvertices(whichdonotexistindefaultRAMSES),making
wheretheα ’saregivenby itthereforemorecomputationallyefficient. InSec.VID,we
i
compareresultsbasedonEqs.(13)and(18).
α =Q ,
1 1
α =Q −Q ,
2 2 1
α =Q −Q , H. Fieldvaluesatcellvertices
3 3 1
α =Q −Q ,
4 5 1
In default RAMSES, the field values are evaluated at cell
α =Q −Q −Q +Q ,
5 4 3 2 1 centers, but the application of Eq. (18) above requires them
α6 =Q7−Q5−Q3+Q1, to be interpolated to cell vertices. When doing so, one can
α =Q −Q −Q +Q , ensurethatthefieldsreconstructedwithtrilinearinterpolation
7 6 5 2 1
(cf.Eq.(14))varycontinuouslyfromcelltocellandalsothat
α =Q −Q −Q −Q +Q +Q +Q −Q ,(15)
8 8 7 6 4 2 5 3 1 thetotalmass11ispreserved. Thisistrivialonaregularmesh
10ThepointsχAandχB canlieanywhereinsidethecellandnotnecessar-
ilyattheintersectionoftheline-of-sightwiththecellfaces. Thisis,for 11Inthissubsection,weusetheword"mass"todenotetheproductofasur-
instance,thegeneralcaseatthestartandendofparticletimesteps. faceorvolumedensitywithanareaorvolume.
10
theyaresharedbybothcoarseandfinecells. Letusconsider
vertex V3 as an example, which is shared by cells c1, c4, f2
andf5. ItisnaturaltoassumethatthedensityatV3depends
onthedensityineachofthesefourcells. Thesimplestwayto
writethisisas
ρ =α ρ +α ρ +α ρ +α ρ , (19)
V3 c1 c1 c4 c4 f2 f2 f5 f5
wheretheα’saresomeweightstobedetermined.Letusfocus
onα . Themassassociatedwithcellc1isρ h2,wherehin
c1 c1
thissubsectionisthecoarsecellsize. Aftertheinterpolation,
wewantthemassofthiscelltobethesumofthemassesas-
sociatedwitheachofitsvertices. Hence,themassassociated
withvertexV3duetocellc1ismc1 =ρ h2/4. Duetomass
V3 c1
conservation,thevalueofmc1 isre-distributedbyvertexV3
V3
toothercellsattheboundary. Inthissense,wecancolloqui-
allydescribevertexV3asamassreservoirthatiscollecting
massfromc1andredistributingittoneighbouringcells. This
massdistributionconstraintcanbewrittenas
ρ h2 α ρ (cid:32) (cid:18)h(cid:19)2 (cid:18)h(cid:19)2
c1 = c1 c1 h2+h2+ +
4 4 2 2
FIG. 6. Two-dimensional illustration of a boundary region. The 1(cid:18)h(cid:19)2 1(cid:18)h(cid:19)2 1(cid:18)h(cid:19)2 1(cid:18)h(cid:19)2(cid:33)
coarse cells are labelled from c1, c2, etc., and a few fine cells are + + + + .
labelled f1, f2, etc. The vertices V1 to V6 are illustrative vertices 2 2 2 2 2 2 2 2
wherethefieldsvaluesaretobeevaluatedfrominterpolationfrom
(20)
thecellcenters,whichiswhereRAMSESevaluatesthefieldsbyde-
fault. Inthisequation,thefirsttwotermsontheRHSrepresentthe
mass from V3 due to cell c1 that is redistributed to cells c1
andc4(α ρ /4isthedensitythatV3contributestoc1and
c1 c1
without AMR, but requires special care at the boundary of c4andh2 istheareaofcellsc1andc4). Thethirdandforth
refinedregions. Wediscussthesesubtleissuesinthissubsec- termsarethesameasthefirsttwo,butforcellsf2andf5(note
tion. that (h/2)2 is the area of the fine cells). The last four terms
Figure 6 illustrates a refinement boundary in two dimen- on the RHS of Eq. (20) must be included as vertex V3 also
sions(weshalllatergeneralizethediscussionstothreedimen- contributestothemassesincellsf1,f2,f5andf8viavertices
sions). Inthefigure,theninecoarselevelcellsshownarela- V2 and V6. Since the contribution from V3 to ρV2 and ρV6
belledfromc1toc9andafewfinelevelcellsarelabelledfrom is only ρV3/2, each of these four terms gets a factor of 1/2
f1 to f12. The vertices V1 to V6 represent different types of comparedtothethirdandforthterms.Wecansolveforαc1in
verticesatwhichonewishestheevaluatethefields.Tobecon- Eq.(20)and,tofacilitatethediscussionsbelow,wewritethe
crete,weshalltaketheexampleofthematterdensityfield,ρ, resultas
buttheinterpolationschemethatwedescribehereisapplied  −1
toanytypeoffield. 1(cid:88)Nf
For vertices of the type of vertex V1, which is shared αc1 =Nc+ 4 Nfv,i , (21)
i
by four cells of the same level, one can simply set ρ =
V1
(cid:0) (cid:1)
ρ +ρ +ρ +ρ /4,whereρ arethevaluesattherelevant whereN isthenumberofcoarsecellsthatsharevertexV3,
f5 f6 f8 f9 i c
cell centers. This average ensures that ρ varies continuously N isthenumberoffinecellsthatsharevertexV3andN is
f fv,i
fromonecelltoanotherwhencomputedwithEq.(14)(orits thenumberofverticesthatthefinecellithatsharesV3hasat
two-dimensionalversion).Inregularmeshsimulations,thisis theboundary(i=f2,f5). Explicitly,forV3wehaveN =2,
c
theonlytypeofvertex. ThevertexV2,whichliesatthemid- N = 2 and N = N = 2. The coefficient α in
f fv,f2 fv,f5 c4
dleofacoarsecelledgeandissharedbytwofinecells,canbe Eq. (19) is obtained in the same way: Eq. (20) remains the
straightforwardlycomputedafterthevaluesinverticesV3and same, just with ρ replaced by ρ (which appears on both
(cid:0) (cid:1) c1 c4
V4havebeendetermined. Explicitly,ρ = ρ +ρ /2. sides of the equation and therefore cancels). Hence, α =
V2 V3 V4 c1
This ensures, for instance, that two rays that are infinitesi- α .
c4
mally close to vertex V2, with one crossing cell c4 and an- The determination of the remaining coefficients, α and
f2
othercrossingfinecellsf5orf8,experiencethesamedensity α differs from α and α . To determine α , the RHS
f5 c1 c4 f2
field (as reconstructed with trilinear interpolation), i.e., there ofEq.(20)remainsthesame, justwithρ andα replaced
c1 c1
isnosharpdiscontinuitybetweenthedensityexperiencedby by ρ and α , respectively. From the reasoning that led to
f2 f2
onerayandtheother. Eq.20,onecouldnaivelythinkthattheLHSwouldbesimply
ThecasesofverticesV3,V4andV5aremoreinvolvedas given by the mass that vertex V3 collects from f2, mf2 =
V3