the agora high-resolution galaxy simulations comparison project. ii. isolated disk test

34 Pages·2016·13.54 MB·English

by  Ji-hoon Kim

TheAstrophysicalJournal,833:202(34pp),2016December20 doi:10.3847/1538-4357/833/2/202
©2016.TheAmericanAstronomicalSociety.Allrightsreserved.
THE AGORA HIGH-RESOLUTION GALAXY SIMULATIONS COMPARISON PROJECT. II.
ISOLATED DISK TEST
Ji-hoon Kim1,2,3,36, Oscar Agertz4,5, Romain Teyssier6, Michael J. Butler7,37,38, Daniel Ceverino8,37,38,
Jun-Hwan Choi9,37,38, Robert Feldmann6,10,37, Ben W. Keller11,37,38, Alessandro Lupi12,37,38, Thomas Quinn13,37,38,
Yves Revaz14,37,38, Spencer Wallace13,37, Nickolay Y. Gnedin15,16,17,38, Samuel N. Leitner18,38, Sijing Shen19,38,
Britton D. Smith20,38, Robert Thompson21,38, Matthew J. Turk22,38, Tom Abel1,2, Kenza S. Arraki23,
Samantha M. Benincasa11, Sukanya Chakrabarti24, Colin DeGraf19, Avishai Dekel25, Nathan J. Goldbaum21,
Philip F. Hopkins3, Cameron B. Hummels3, Anatoly Klypin23, Hui Li26, Piero Madau12,27, Nir Mandelker25,28,
LucioMayer6,KentaroNagamine29,30,SarahNickerson6,BrianW.O’Shea31,JoelR.Primack32,SantiRoca-Fàbrega25,
VadimSemenov16,IkkohShimizu29,ChristineM.Simpson33,KeitaTodoroki34,JamesW.Wadsley11,andJohnH.Wise35
(for the AGORA Collaboration)39
1KavliInstituteforParticleAstrophysicsandCosmology,SLACNationalAcceleratorLaboratory,MenloPark,CA94025,USA;[email protected]
2DepartmentofPhysics,StanfordUniversity,Stanford,CA94305,USA
3DepartmentofAstronomy,CaliforniaInstituteofTechnology,Pasadena,CA91125,USA
4DepartmentofPhysics,UniversityofSurrey,Guildford,Surrey,GU27XH,UK
5LundObservatory,DepartmentofAstronomyandTheoreticalPhysics,LundUniversity,SE-22100Lund,Sweden
6CentreforTheoreticalAstrophysicsandCosmology,InstituteforComputationalScience,UniversityofZurich,Zurich,8057,Switzerland
7Max-Planck-InstitutfürAstronomie,D-69117Heidelberg,Germany
8ZentrumfürAstronomiederUniversitätHeidelberg,InstitutfürTheoretischeAstrophysik,D-69120Heidelberg,Germany
9DepartmentofAstronomy,UniversityofTexas,Austin,TX78712,USA
10DepartmentofAstronomy,UniversityofCaliforniaatBerkeley,Berkeley,CA94720,USA
11DepartmentofPhysicsandAstronomy,McMasterUniversity,Hamilton,ONL8S4M1,Canada
12Institutd’AstrophysiquedeParis,SorbonneUniversites,UPMCUnivParis6etCNRS,F-75014Paris,France
13DepartmentofAstronomy,UniversityofWashington,Seattle,WA98195,USA
14InstituteofPhysics,Laboratoired’Astrophysique,ÉcolePolytechniqueFédéraledeLausanne,CH-1015Lausanne,Switzerland
15ParticleAstrophysicsCenter,FermiNationalAcceleratorLaboratory,Batavia,IL60510,USA
16DepartmentofAstronomyandAstrophysics,UniversityofChicago,Chicago,IL60637,USA
17KavliInstituteforCosmologicalPhysics,UniversityofChicago,Chicago,IL60637,USA
18DepartmentofAstronomy,UniversityofMaryland,CollegePark,MD20742,USA
19KavliInstituteforCosmology,UniversityofCambridge,Cambridge,CB30HA,UK
20InstituteforAstronomy,UniversityofEdinburgh,RoyalObservatory,EdinburghEH93HJ,UK
21NationalCenterforSupercomputingApplications,UniversityofIllinois,Urbana,IL61801,USA
22SchoolofInformationSciences,DepartmentofAstronomy,UniversityofIllinois,Urbana,IL61801,USA
23DepartmentofAstronomy,NewMexicoStateUniversity,LasCruces,NM88001,USA
24SchoolofPhysicsandAstronomy,RochesterInstituteofTechnology,Rochester,NY14623,USA
25CenterforAstrophysicsandPlanetaryScience,RacahInstituteofPhysics,TheHebrewUniversity,Jerusalem,91904,Israel
26DepartmentofAstronomy,UniversityofMichigan,AnnArbor,MI48109,USA
27DepartmentofAstronomyandAstrophysics,UniversityofCaliforniaatSantaCruz,SantaCruz,CA95064,USA
28DepartmentofAstronomy,YaleUniversity,NewHaven,CT06520,USA
29DepartmentofEarthandSpaceScience,GraduateSchoolofScience,OsakaUniversity,Toyonaka,Osaka,560-0043,Japan
30DepartmentofPhysicsandAstronomy,UniversityofNevada,LasVegas,NV89154,USA
31DepartmentofComputationalMathematics,ScienceandEngineering,DepartmentofPhysicsandAstronomy,NationalSuperconducting
CyclotronLaboratory,MichiganStateUniversity,Lansing,MI48824,USA
32DepartmentofPhysics,UniversityofCaliforniaatSantaCruz,SantaCruz,CA95064,USA
33HeidelbergerInstitutfürTheoretischeStudien,D-69118Heidelberg,Germany
34DepartmentofPhysicsandAstronomy,UniversityofKansas,Lawrence,KS66045,USA
35CenterforRelativisticAstrophysics,SchoolofPhysics,GeorgiaInstituteofTechnology,Atlanta,GA30332,USA
Received2016October7;revised2016October25;accepted2016October28;published2016December19
ABSTRACT
Using an isolated Milky Way-mass galaxy simulation, we compare results from nine state-of-the-art gravito-
hydrodynamicscodeswidelyusedinthenumericalcommunity.Weutilizetheinfrastructurewehavebuiltforthe
AGORA High-resolution Galaxy Simulations Comparison Project. This includes the common disk initial
conditions, common physics models (e.g., radiative cooling and UV background by the standardized package
GRACKLE) and common analysis toolkit yt, all of which are publicly available. Subgrid physics models such as
Jeans pressure floor, star formation, supernova feedback energy, and metal production are carefully constrained
across codeplatforms. With numerical accuracy that resolves thedisk scale height,we find that thecodes overall
agree well with one another in many dimensions including: gas and stellar surface densities, rotation curves,
36 EinsteinFellow.
37 Theseauthors,inalphabeticalorder,contributedtothearticlebyleadingthe
effortwithineachcodegrouptoperformandanalyzesimulations.
38 Theseauthors,inalphabeticalorder,contributedtothearticlebydeveloping
GRACKLEandimplementingitsinterfacewithparticipatingcodes.
39 http://www.AGORAsimulations.org/
1
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
velocity dispersions, density and temperature distribution functions, disk vertical heights, stellar clumps, star
formation rates, and Kennicutt–Schmidt relations. Quantities such as velocity dispersions are very robust
(agreementwithinafewtensofpercentatallradii)whilemeasureslikenewlyformedstellarclumpmassfunctions
show more significant variation (difference by up to a factor of ∼3). Systematic differences exist, for example,
between mesh-based and particle-based codes in the low-density region, and between more diffusive and less
diffusive schemes in the high-density tail of the density distribution. Yet intrinsic code differences are generally
smallcomparedtothevariationsinnumericalimplementationsofthecommonsubgridphysicssuchassupernova
feedback. Our experiment reassures that, if adequately designed in accordance with our proposed common
parameters, results of a modern high-resolution galaxy formation simulation are more sensitive to input physics
than to intrinsic differences in numerical schemes.
Keywords:cosmology:theory–galaxies:evolution–galaxies:formation –galaxies:kinematicsanddynamics–
ISM: structure – methods: numerical
1. INTRODUCTION Collaboration has aimed to compare galaxy-scale numerical
experimentsonavarietyofcodeplatformswithstate-of-the-art
Decadesofstrenuouseffortbycomputationalastrophysicists
resolution. Our shared goal is to ensure that physical
have propelled numerical experiments to become one of the
assumptions are responsible for any success in galaxy
most widely used tools in theorizing how galaxies form in the
formation simulations, rather than artifacts of particular
universe. Numerical experiments are often the only means to
implementations. Through a multi-platform approach from
put our theory to a test, the result of which we can compare
the beginning, we strive to improve all our codes by
with observational data to validate the model’s feasibility. “increasingthelevelofrealismandpredictivepowerofgalaxy
Since the success of galaxy formation theory is predicated on
simulations and the understanding of the feedback processes
robust numerical experiments, it is only reasonable that we that regulate galaxy metabolism” (Kim et al. 2014), and by
apply the same scientific standard of reproducibility to galaxy doing so to find solutions to long-standing problems in galaxy
formation simulations. In other words, it should be considered
formation. Because the interplay between numerical resolution
as a fundamental principle that researchers must not establish
and subgrid modelings of stellar physics is crucial in galaxy-
findings from a single numerical experiment as scientific
scalesimulations,werequirethatsimulationsbedesignedwith
knowledge. Only after the result is reproduced independently state-of-the-artresolution,100pc,whichiscurrentlyallowed
byotherresearchersandprovennottobeanisolatedincidence
within realistic computational cost bounds.
canwebuildanyconclusivetheoryabouthowgalaxiesactually In the Project’s flagship paper (Kim et al. 2014), we
form in the universe. explained the philosophy behind the Project and detailed the
However, the task of replicating galaxy simulations or, publicly available Project infrastructure we have put together.
equivalently, comparing simulations between codes, has not Wealsodescribedtheproof-of-concepttest,inwhichwefield-
received high priority.40 Instead, the task is considered tested our infrastructure with a dark matter-only cosmological
complex and time-consuming because one needs to ensure zoom-in simulation, finding a robust convergence between
that identical physics is used in an identical initial condition participating codes. More than 140 researchers from over 60
(IC) with identical runtime settings. This is sometimes academicinstitutionsworldwidehavesinceagreedtotakepart
perceived as tedious and unrewarding for early-career in the Collaboration, many of whom having been actively
researchers. In fact, the lack of reproducibility checks is not engaged in working groups and sub-projects.41 The cohort of
unique to the field of numerical galaxy formation (e.g., Open numerical codes participating in the Project currently include,
Science Collaboration 2015; Nature Survey 2016). And its but are not limited to in future studies: the Lagrangian
cause is not simply an unwillingness of only computational smoothed particle hydrodynamics codes (SPH; Gingold &
astrophysicists, either (Everett & Earp 2015). Rather, Monaghan 1977; Lucy 1977; Monaghan 1992) CHANGA,
addressing the system (or the lack thereof) which checks the GADGET, GASOLINE, and GEAR, and the Eulerian adaptive
reproducibilityofsimulationswouldrequireacollectiveaction mesh refinement codes (AMR; Berger & Oliger 1984; Berger
by the entire community. It cannot be simply about asking & Colella 1989) ART-I, ART-II, ENZO, and RAMSES, and the
individual researchers to release their data dumps, but should mesh-free finite-volume Godunov code GIZMO (see Section 5
be about building a system that incentivizes simulations for information on each code).
published in a reproducible manner. It should also be about In this second report of our continuing endeavor, we use
assembling an infrastructure that reduces the cost of reprodu- an isolated Milky Way-mass galaxy simulation to compare
cibilitychecks,onwhichsimulationsareverifiedroutinelyand ninewidelyusedstate-of-the-artgravito-hydrodynamicscodes.
effortlessly (Begley & Ioannidis 2015; Nosek et al. 2015). As in all comparison studies in AGORA, the participating
The AGORA High-resolution Galaxy Simulations Compar- codes share a common IC (i.e., generated by MAKEDISK;
isonProject(AssemblingGalaxiesOfResolvedAnatomy)isthe see Section 2), common physics models (e.g., radiative
collective response by the numerical galaxy formation com- cooling and UV background provided by the standardized
munity to such a challenge. Since its first meeting in 2012 at package GRACKLE; see Section 3.1; Bryan et al. 2014; Kim
the University of California at Santa Cruz, the AGORA etal.2014;Smithetal.2016),42andcommonanalysisplatform
40 Code comparisons in the astrophysical community have previously been 41 SeetheProjectwebsiteathttp://www.AGORAsimulations.org/formore
undertaken, albeit with simplified physics in a different scale (e.g., Frenk informationontheProjectincludingitsmembership,anditstask-orientedand
etal.1999;O’Sheaetal.2005),orfocusingonlyonhydrodynamicssolvers science-orientedworkinggroups.
(e.g.,Agertzetal.2007;Taskeretal.2008). 42 Thewebsiteishttp://grackle.readthedocs.org/.
2
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
Table1
InitialConditionCharacteristicsa
DarkMatterHalo StellarDisk GasDisk StellarBulge
Densityprofile Navarroetal.(1997) Exponential Exponential Hernquist(1990)
Structuralproperties M200=1.074´1012M,vc,200=150kms-1, Md,=3.438´1010M, Md,gas=8.593´109M, Mb,=4.297´109M,
R200=205.5kpc,c=10,l=0.04 rd=3.432kpc,zd=0.1rd fgas =0.2 Mb, Md=0.1
Numberofparticles 105 105 105 1.25´104
Particlemass mDM=1.254´107M m,IC=3.437´105M mgas,IC=8.593´104M m,IC=3.437´105M
Note.
aFordetailedexplanationsontheseparameters,seeSection2.
Figure1.0MyrsnapshotsoftheisolatedMilkyWay-massgalaxysimulationsbynineparticipatingcodes.Diskgassurfacedensitiesina30kpcbox,edge-on(top)
andface-on(bottom),producedwiththecommonanalysistoolkityt.Forvisualizationsoftheparticle-basedcodeshereafter(Figures1–3,14–15,32,34,and35)—
butnotinanyotheranalysesexceptthesefigures—ytusesanin-memoryoctreeonwhichgasparticlesaredepositedusingsmoothingkernels.Comparing0Myr
snapshots—dumpedimmediatelyaftereachcodereadsintheIC—istochecktheexactidentityofICsinterpretedbyeachcode.SeeSection6.1formoreinformation
onthisfigure,andSection5fordescriptionsofparticipatingcodesinthiscomparison.Thehigh-resolutionversionsofthisfigureandarticleareavailableattheProject
website,http://www.AGORAsimulations.org/.
(i.e.,yttoolkit;Turketal.2011).43Weadoptspatialresolution in the numerical galaxy formation community the AGORA
of 80 pc that resolves the scale height of the disk. This helps Project strives to promote.
the codes to be less dependent on phenomenological prescrip-
tions of sub-resolution processes which are inevitably intro- 2. INITIAL CONDITION
duced in low-resolution (>kpc) simulations. As modern
In this section we describe the Milky Way-mass isolated IC
galaxy formation simulations with state-of-the-art resolution weadoptinthisstudy.WhilethisICispartofasetofdiskICs
and physics prescriptions become more and more computa- generated for AGORA simulations that were first introduced in
tionallyexpensive,itistimelythatwecomparehigh-resolution Section 2.2of the Project flagship paper (Kim et al.2014), we
isolated disk simulations to check how successfully these briefly explain its important structural properties for
galaxies are reproduced by their peers.44 Readers should note completeness.45
that our intention is not to identify a “correct” or “incorrect” The disk galaxy IC with properties characteristic of Milky
code,buttofocusinsteadonjuxtaposingthecodesforphysical Way-mass galaxies at redshift z ~ 1 is generated with a
insights and learn how much scatter one should expect among privately shared version of MAKEDISK (Springel et al.
modern numerical tools in the field (see Section 7 for more 2005).46,47 The IC has the following components (see also
discussion). Table 1 and Figure 1): (1) a dark matter halo with
The remainder of this paper is organized as follows. In M200 = 1.074 ´ 1012M, R200 = 205.5kpc and circular velo-
Section2weexplaintheisolateddiskICusedinthestudy.The city of vc,200 = 150kms-1 that follows the Navarro–Frenk–
White (Navarro et al. 1997) profile with concentration
common input physics and runtime parameters required in the
parameter c=10 and spin parameter l = 0.04; (2) an
participating codes are discussed in Sections 3 and 4,
exponential disk with M = 4.297 ´ 1010 M , scale length
respectively. Then Section 5 describes nine hydrodynamics d 
r = 3.432kpcandscaleheightz = 0.1r thatiscomposedof
codes that participated in this comparison. Section 6 presents d d d
80%starsand20%gasinmass(i.e., f = M M = 0.2);
the results of our comparison focusing on similarities and gas d,gas d
(3) a stellar bulge with M = 4.297 ´ 109 M that follows
discrepancies discovered in various multi-dimensional ana- b, 
lyses. Finally in Section 7 we summarize our findings and the Hernquist profile (Mb, Md = 0.1; Hernquist 1990). The
conclude the paper with remarks on future work. We will also
45 ThepublicDropboxlinkishttp://goo.gl/8JzbIJ.
stresstheimportanceofcollaborativeandreproducibleresearch
46 MAKEDISK is an earlier realization of a code similar to GALIC (Yurin &
Springel2014).GALICispubliclyavailable,anditswebsiteishttp://www.h-
43 Thewebsiteishttp://yt-project.org/. its.org/tap-software-en/galic-code/.
44 Comparisonsofcosmologicalzoom-insimulationsarealsointhemakingto 47 WhiletheMilkyWay’s f is∼10%,typicalgalaxieswiththeMilkyWay
test the robustness of the code suite over 13.8Gyr of evolution. See the stellarmassatz~0 have fgas ~20%(Catinellaetal.2010).Inthisregard,
Project’sflagshippaper(Kimetal.2014)formoreinformation. onecansaythatwemodelamgasoretypicalgalaxyatz~0thantheMilkyWay.
3
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
disk or bulge stars in the IC do not contribute to the feedback cooling, star formation and feedback. In Section 3.1, we first
budget. listthegasphysicsthatarecommoninbothSim-noSFandSim-
Among the three resolution choices provided in Kim et al. SFF. Then in Section 3.2, the subgrid prescriptions of stellar
(2014), here we employ a “low-resolution” IC which has 105 physics for Sim-SFF are explained.
particleseachforthehalo,thestellardisk,andthegasdisk,and
1.25 ´ 104 particles for the bulge. The initial gas temperature 3.1. Gas Physics: Radiative Cooling, UV Background, and
in the disk is set to 104 K, not to the specific internal energy Pressure Floor
computed by MAKEDISK. The initial metal fraction in the gas
The rate at which the gas in our galaxy radiatively cools is
disk is 0.02041.48 When the gas disk is initialized on mesh-
determined by AGORAʼs standard chemistry and cooling
based codes (ART-I, ART-II, ENZO, and RAMSES), instead of library GRACKLE (Bryan et al. 2014; Kim et al. 2014; Smith
usingtheparticlesprovidedbyMAKEDISK,werequirethatthe et al. 2016, see footnote 42). For this study, the equilibrium
participants use an analytic density profile of
coolingversionofGRACKLEisinterfacedwitheachparticipat-
r (r, z) = r e-r rd · e-∣z∣ zd (1) ing code, either via GRACKLEʼs original interface or via N.
d,gas 0 Gnedin’s auxiliary API.50 In the chosen equilibrium cooling
with r = M (4pr2z ), where r is the cylindrical radius mode,GRACKLE follows tabulated coolingratespre-computed
0 d,gas d d by the photoionization code CLOUDY (Ferland et al. 2013).51
andzistheverticalheightfromthediskplane.Tosetupadisk
The pre-computed look-up table also includes metal cooling
in a centrifugal equilibrium we also ask that the participants
rates for solar abundances,1Z , as a function of gas number
utilize the rotational velocity profile binned from an actual gas 
density and temperature. These metal cooling rates are then
particledistributionwithin∣z∣ < zd.49Inmesh-basedcodes,we scaled linearly with metallicity which is followed in our
additionally include a uniformly low-density gas halo with simulations as a separate passive scalar.52 We also adopt
n = 10-6cm-3 and zero initial velocity, since they cannot metagalactic UV background radiation at z=0 by Haardt &
H
have cells with zero density. The halo is initially set to 106 K Madau (2012) provided by GRACKLE. For the difference
andzerometallicity. Notethatthisgaseoushalodoesnotexist between the chosen UV background model and previous
in particle-based codes (SPH codes or GIZMO). calculations such as Haardt & Madau (1996) or Faucher-
It is worth noting one point about the common disk IC Giguèreetal.(2009),wereferthereaderstoSection3.3ofKim
adopted here. As readers may find in Section 6.6, in our et al. (2014).
experiment that includes radiative cooling, star formation and Lastly, a non-thermal Jeans pressure floor is applied to
feedback(“Sim-SFF”;seeSection3),only~109M additional stabilize the scales of the smoothing length (particle-based
starsformin500Myrinallcodesonaverage.Whencompared codes)orthefinestcell(mesh-basedcodes)againstunphysical
with the stellar and gas components present in the IC, this collapseandtoavoidartificialfragmentationduetounresolved
means only a ~3% increase in stellar mass, and a ~12% pressure gradient (Truelove et al. 1997; Robertson &
decrease in gas mass. As the discussion in Section 6 should Kravtsov 2008). In practice, it is achieved by enforcing that
makeclear,giventhisrelativelysmallchangeinstellarmass,it the local Jeans length lJeans be sufficiently resolved with the
is expected that the stellar feedback in our experiment will be finest resolution elements at all times. That is,
very inefficient.
l = N Dx (2)
Jeans Jeans
3. COMMON PHYSICS: Sim-noSF AND Sim-SFF whereDx = 80pc is the adopted spatial resolution (finest cell
size or softening length; see Section 4.1) and N = 4 is the
We now describe the common physics employed in our Jeans
Jeans number adopted from Truelove et al. (1997). This gives
experimentincludingequilibriumgascooling,metagalacticUV
background, star formation, and energy and metal yields by the required pressure floor value as
supernovae. Note that the common physics adopted here is a
1
variation of the common physics model recommended in all P = N2 Gr2 Dx2 (3)
Jeans gp Jeans gas
AGORA simulations by default; see Section 3 of Kim et al.
(2014). For the present study all participating code groups are
where G is the gravitational constant,g = 5 3is the adiabatic
askedtoruntwosimulationsstartingfromtheidenticalIC:(1)
index, and r is the gas density. Note that N is not
“Sim-noSF” with radiative gas cooling but without star gas Jeans
formation or feedback, and (2) “Sim-SFF” with radiative necessarily equal to the parameter controlling the pressure
supportineachcode.Foractualparameterchoicesforselected
48 Thisfractionalvalue0.02041correspondsto1ZforGRACKLEv2.0,butto codes, see Appendix A. For implementations using polytropes
1.5761 ZforGRACKLEv2.1orabove.Itisbecausethesolarmetallicityunit in ART-II and RAMSES, see Sections 5.2 and 5.4, respectively.
Z was updated from 0.02041 to 0.01295 in GRACKLE v2.1. Since cooling
rthateescporoel-itnagbulraatteedsbyinCLGORUADCYKaLrEeʼsat1eqZuil,ibnroitumatspceocoilfiincgmemtaoldferacwtiiollnvdailfufeer, 3.2. Stellar Physics: Star Formation, and Energy, Mass, and
depending on which GRACKLE version is adopted (see Section 3.1 for more Metal Yields from Core-collapse Supernovae
onGRACKLE).Forexample, inthe current study,the codesusingGRACKLE
v2.1 (CHANGA,GASOLINE, GADGET-3, and GIZMO) show slightly enhanced In addition to the gas physics described in the previous
cooling rates than the ones using GRACKLE v2.0 or below (ART-I, ART-II, section, Sim-SFF incorporates subgrid models for star
ENZO,RAMSES,andGEAR).Generallyspeaking,initialgasmetallicityshould
wbeitshetthuepcsoodet.haWteitriesfecroninstiesrteesnttewdirtehadtheerscthootsheenGGRRAACCKKLLEEvv2e.r1siroenleiansteernfoacteinagt 50 Thewebsiteishttps://bitbucket.org/gnedin/agora_api/.
https://goo.gl/BNRfwJ. 51 Thewebsiteishttp://www.nublado.org/.
49 This actual initial velocity profile, provided in vcirc_SPH.dat in our 52 See, however, footnote 48 on how a different version of GRACKLE may
public Dropbox link, is different from the file vcirc.dat produced by affectthecoolingratesforthegaswiththesamemetalfraction(butnotthe
MAKEDISKitself.Thedifferenceis~5%inthecentralfewkpc. samemetallicityinterpretedbyCLOUDY).
4
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
formationandsupernovafeedback.First,aparcelofgasabove 4.1. Gravitational Softening Length and Finest Mesh Size
the threshold n = 10cm-3= r m produces stars
H,thres gas,thres H Forallcodesthegasmassresolutioninhydrodynamicsneeds
at a rate that follows the local Schmidt law as
to be set as close as possible to m = 8.593 ´ 104M .
gas,IC 
Assuming that we wish to resolve a self-gravitating clump with
dr = rgas (4) 64 of these resolution elements, the corresponding Jeans length
dt tff scale becomes
where r is the stellar density, t = (3p (32Gr ))1 2 is the ⎛ 64m ⎞1 3
 ff gas l = 2⎜ gas,IC ⎟ = 348.7pc, (5)
local free-fall time, and = 1% is thestar formation efficiency Jeans ⎝(4p 3)nH,thres⎠
per free-fall time. We caution that the parameter  is not

necessarilyequaltothestarformationefficiencyparameterfound andtherefore,fromEquation(2)wechooseaspatialresolutionof
ineachcode(e.g.,forCHANGAandGASOLINE,seeAppendixB). 80pc.ThisvalueisusedasthefinestcellsizeDx formesh-based
Foranewstarparticletospawn,itshouldhaveatleastthemass codes, and as the gravitational softening length for particle-
grav
ofagasparticleintheIC,m = 8.593 ´ 104M .Notethat basedcodes.Forallparticle-basedcodestakingpartinthepresent
gas,IC 
n adopted in this experiment is for this particularrunonly, study,gravityissoftenedaccordingtothecubicsplinekernel(e.g.,
H,thres
andrepresentswheretheJeanspolytropeintersectswithatypical Equation(A1)ofHernquist&Katz 1989). Forreaders interested
T - r equationofstateinourdisks(seeSections5.2and5.4).53 in the actual parameter choices, in Appendix C we examine the
New star particles inject energy, mass, and metals back into meaningsofrelevantparametersindifferentparticle-basedcodes.
the interstellar medium (ISM) through core-collapse (Type II)
supernovae. Assuming the AGORA standard Chabrier (2003) 4.2. Minimum Hydrodynamical Smoothing Length
initial massfunction(IMF)andthatstarswithmassesbetween For particle-based codes (including GIZMO; see Section 5.9
8 and 40 M explode as Type II supernovae, one Type II and footnote 67), we require that the hydrodynamical
supernova occurs per every 91M stellar mass formed (see
 smoothing lengths for collisional particles do not drop below
Section 3.5 of Kim et al. 2014). With the AGORA
20% of the gravitational softening lengths. Unlike the
recommended fitting formulae Equations (4)–(6) of Kim
gravitational softening kernel, exact smoothing kernel choices
et al. (2014) and the assumed IMF, this single burst is found
differ from code to code, and are detailed for each of the
to release 2.63M of metals54 and14.8M of gas (including
  particle-based codes in Section 5. We also refer the readers
metals).Perevery91M stellarmass,thesemetalandmassare
 interestedintheactualparameterchoicestoAppendixCagain.
instantaneously deposited into its surrounding after a delay
time of 5Myr, along with a net thermal energy of 1051erg. 4.3. Refinement Strategy
We note that exact deposit schemes for energy, mass, and
metalsareleftateachparticipant’sdiscretion.Wedonotintend We recommend to mesh-based code groups that a cell be
to overly specify a single common deposit scheme which will split into 8 child cells once the cell contains more mass than
need to be inevitably different from one code to another (e.g., mgas,IC = 8.593 ´ 104M (1 gas particle mass in the IC of
between mesh-based codes and particle-based codes), as we particle-based codes), or 8 collisionless particles (disk/bulge
arguedinSection3.8ofKimetal.(2014).Nevertheless,forall star particles in the IC with m,IC = 3.437 ´ 105M, or dark
mesh-based codes (ART-I, ART-II, ENZO, and RAMSES), the matterparticleswithmDM = 1.254 ´ 107M).Thiscausesthe
same strategy was chosen: thermal energy, mass, and metals grids to be refined in a fashion similar to the Lagrangian
behavior of particle-based codes, and keeps the ratio of
are added to the cell where a 5Myr old star particle sits at the
collisionless particle numbers to gas cells approximately unity
time of explosion, and to this cell only. For particle-based
codes (CHANGA, GASOLINE, GADGET-3, GEAR, and GIZMO), onaverage.However,exactrefinementstrategiesdifferslightly
each code’s deposit scheme is discussed in detail in Section 5. fromcodetocode,andaredetailedforeachofthemesh-based
codes in Section 5.55 We continue to refine the grids down to
In future AGORA projects, we plan to calibrate different
the resolution limit Dx = 80 pc (see Section 4.1) where the
feedback schemes against observations and against one
non-thermal pressure floor kicks in (see Section 3.1).
another. We refer the readers to Section 7 for more discussion
on this future work.
5. PARTICIPATING CODES
Inthissection we introducetheninegravito-hydrodynamics
4. COMMON RUNTIME PARAMETERS codes taking part in this test, focusing in particular on
hydrodynamics solvers, refinement schemes for mesh-based
Here we review the runtime parameters each group is
codes (ART-I, ART-II, ENZO, and RAMSES), and supernova
required to adopt, such as gravitational softening and
feedback implementations for particle-based codes (CHANGA,
hreyfidnreomdyennatmthicressmhoolodtshfinogr mleensght-hbsasfoedr pcaordtiecsl.e-based codes and GASOLINE, GADGET-3, GEAR, and GIZMO). We leave out
detailsthatarecommonlyadoptedacrossplatformssuchasgas
cooling (Section 3.1) or star formation (Section 3.2), or that
53 AsnotedinKimetal.(2014),starformationprescriptionparameterssuchas were included in the AGORA flagship paper such as
n or, the initial mass of star particles, and the stochasticity of star
H,thres 
formation,areallhighlydependentonnumericalresolution.Anidealizedtest
likethedisksimulationpresentedhereisessentialtotuneupsuchparameters 55 Givendifferencesinrefinementmachineriesamongmesh-basedcodesitis
forcomputationallyexpensivecosmologicalsimulations. impractical, if not impossible, to impose an exactly identical refinement
54 Perunitstellarmassformed,thetotalfractionalejectedmetalmasses(oxygen criterionacrossallcodes.Weinsteadadoptatrial-and-error approachwithin
andironcombined)is MZ =2.09MO+1.06MFe=2.09´0.0133+1.06´ the guideline presented in Section 4.3, which resulted in all mesh-codes
0.0011=2.9%. eventuallyconvergingtoasimilaroverallgridstructure.
5
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
gravitationaldynamics(Section5ofKimetal.2014).Wealso where P is the value entering the Riemann solver, P is the
cell gas
point out that the codes involved in future AGORA studies are gaspressurefieldinthesimulation,k istheBoltzmannconstant,
B
not necessarily limited to the ones described herein. n = r m , andT = T (n n ) with T = 1800K and
H gas H Jeans J H H,J J
n = 8cm-3. This polytrope choice is designed to match the
5.1. ART-I H,J
common prescription Equation (3) with N  4. For the
Jeans
In ART-I, differential equations of fluid dynamics are supernova feedback scheme to deposit the thermal energy
integrated using a shock-capturing Eulerian method described
adopted by all mesh-based codes, see Section 3.2.
in Khokhlov (1998). It uses a second-order accurate Godunov
solver (Godunov 1959) that evaluates Eulerian fluxes by
solvingtheRiemannproblemateverycellinterface(Colella& 5.3. ENZO
Glaz 1985). Left and right states of the Riemann problem are ENZO is a block-structured adaptive mesh code, developed by
obtained by piecewise linear interpolation (van Leer 1979). In an open-source, community-driven approach (Bryan & Nor-
contrast to other versions of ART (see Section 5.2.1 of Kim man 1997; O’Shea et al. 2004; Bryan et al. 2014).57 Among a
etal.2014),ART-Iwithdistinctivestarformationandfeedback variety of solver choices, for this comparison the third-order
recipes (e.g., Ceverino & Klypin 2009; Ceverino et al. 2014) accuratepiecewiseparabolicmethodisselectedtoreconstructthe
have been developed by A. Klypin and collaborators. left and right states of the Godunov problem (Colella &
The octree-based, multi-level adaptive mesh allows users to Woodward 1984; Bryan et al. 1995), along with a Harten–Lax–
control the grid structure at the individual cell level. For this vanLeerwithcontact(HLLC)Riemannsolver(Toroetal.1994).
comparison, the ART-I group uses a 1283 root grid covering a A maximum 30% of the required Courant–Friedrichs–Lewy
(1.304Mpc)3 box, then achieves an ∼80 pc cell size at (CFL)timestepisusedtoadvancefluidelements;i.e.,CFLsafety
maximum 7 levels of refinement. The mass thresholds above factor=0.3.Inadditiontosolvingtheconservationequationsfor
whichacellisadaptivelyrefinedintoanoctof8childcellsare mass, momentum and energy, the equation for internal energy is
mgas,IC = 8.593 ´ 104M and m,IC = 3.437 ´ 105M for also solved in parallel, and the conservative or non-conservative
gas and collisionless particles, respectively.56 For the super- formulationisadaptivelyselectedbasedonalocalestimateofthe
nova feedback scheme to deposit the thermal energy adopted energy truncation errors. This ensures that the gas temperature
by all mesh-based codes, we refer the readers to Section 3.2. remains physical, even in highly supersonic regions.
The ENZO group uses a 643 initial root grid covering a
5.2. ART-II (1.311Mpc)3simulationbox,thenachieves80pcresolutionwith
maximum 8 levels of refinement. The mass thresholds above
ART-II solves the gravito-hydrodynamics equations using a
particle-mesh+EulerianAMRapproach.ART-IIfeaturesMPI whichacellisrefinedbyfactorsoftwoineachaxisaremgas,IC =
parallelization for distributed memory machines, flexible time- 8.593 ´ 104M and 8m,IC = 8 ´ 3.437 ´ 105M for gas
stepping hierarchy, and a variety of unique physics modules andcollisionlessparticles,respectively(seefootnote56).Thenon-
(e.g., Gnedin & Kravtsov 2011; Agertz et al. 2013) developed thermal pressure floor Equation (3) is used to modify the gas
by N. Gnedin, A. Kravtsov and collaborators. pressureinsidetheRiemannsolver,butnottoaltertheactualgas
For the present study, starting from a uniform 1283 root grid energy field. For the supernova feedback scheme adopted by all
covering(1.311Mpc)3,cellsarerefinedupto7additionallevels mesh-based codes, we refer the readers to Section 3.2.
toreachthefinestsizeof80pc.Sphericalregionsof4(6,10)root
grid cells radius around the box center are always refined to at 5.4. RAMSES
least 3 (2, 1) additional levels relative to the root grid. The (de-)
refinement procedure consists of three steps. First, cells are RAMSES isanoctree-based adaptive meshcodefeaturingan
marked for refinement if the gas mass in the cell exceeds unsplit second-order accurate Monotone Upstream-centered
0.6m = 0.6 ´ 8.593 ´ 104M ,orifthecellcontainstwo Scheme for Conservation Laws (MUSCL) Godunov scheme
ormograes,IdCarkmatterand/orstarparticlesthatwerepresentinthe for the gaseous component (Teyssier 2002).58 For this
IC.Wethenuseadiffusionsteptoalsomarkneighboringcellsfor comparison, RAMSES group uses a ideal gas equation of state
refinementandthussmooththeshapeoftheregionstoberefined. with g = 5 3, along with the HLLC Riemann solver (Toro
By contrast, cells with gas masses below 0.2m or without et al. 1994) and the MinMod slope limiter (Roe 1986). The
gas,IC
particlesaremarkedforde-refinementprovidedtheyalsosatisfya CFLsafetyfactorforcontrollingthetimestepissetto0.5.The
number of additional constraints. Finally, cells are refined (de- dual energy formalism adopted in ENZO simulations
refined) by splitting them into 8 (by merging 8 child cells). (Section 5.3) is also used in RAMSES runs.
ThepressurefloorimplementedinART-IIaffectscellsatthe For this study, starting from a uniform 1283 root grid
highest level of refinement by modifying the gas pressure covering(320kpc)3, cells are refined up to 5 additional levels
values that enter the Riemann solver (i.e., not the actual toachievean∼80pccellsize.Therefinementprocessworksas
pressure or temperature fields) with follows. First, new refinement is triggered on a cell-by-cell
basisif thebaryonicmass (gas+ newlyformed stars) exceeds
Pcell = max(PJeans, Pgas) (6) mgas,IC = 8.593 ´ 104M, or if the number of dark matter
and/orstarparticlesthatarepresentintheICexceeds8.59We
=max(n k T , P ) (7)
H B Jeans gas
57 Thewebsiteishttp://enzo-project.org/.
56 WenotethatART-IandENZOcannotrefinecellsbyparticlenumbers,but 58 Thewebsiteishttp://www.itp.uzh.ch/~teyssier/Site/RAMSES.html.
onlybyparticlemasses.Bycontrast,inthereportedruns,ART-IIandRAMSES 59 Readers should notice subtle differences here in refinement strategies
refinecellsbyparticlenumbers.Therefinementcriteriaarechosentoensurean between RAMSES and other mesh-based codes. Newly formed stars are
agreement among mesh-based codes in overall grid structures. See also considered as part of the baryonic fluid, so they do not change the particle
footnote55. refinementbasedsolelyoncollisionlessparticlesintheIC.
6
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
then mark additional cells by performing a mesh smoothing 5.6. GASOLINE
operation,expandingtheinitialareabyonecellwidthinevery
GASOLINE is a massively parallel SPH code, first described
direction. When new cells are created or old cells destroyed,
in Wadsley et al. (2004), that has subsequently been updated
density, momentum and internal energy are used as averaging
with modern SPH features. It contains a subgrid model for
and interpolating variables, thereby preventing a grid point
turbulentmixingofmetalsandenergy(e.g.,Shenetal.2010),a
with spurious temperature.
In RAMSES, the gas pressure field includes the non-thermal timestep limiter by Saitoh & Makino (2009) (see Section 5.5),
and a geometric density estimator for SPH force expressions
pressure support term given by a temperature polytropeT =
mT (n¢ n ) with mean molecular weight μ, n¢ = r XJeanms , (see Section 2.4 of Keller et al. 2014, for a latest detailed
J H H,J H gas H H description of the code and its performance).
TJ = 1800K, nH,J = 8cm-3 and XH = 0.76.60 As in ART-II, Forthecurrentwork,theGASOLINEgroupusesaWendland
this polytrope approximately matches the common pressure C4smoothingkernel(Dehnen&Aly2012)withN = 200
support prescription Equation (3). Newly created star particles in neighbors.64 The same feedback scheme assmCooHthANGAʼs
RAMSEShaveafixedmassofmgas,IC,buttheyarespawnedwitha (Section 5.5) is implemented, smoothed with the Wendland
Poissonprobabilitydistributionwhoseparametersaredesignedto C4 kernel over 64 neighbors (not N = 200) to better
mimic the local Schmidt law, Equation (4). Lastly, for the smooth
match the amount of mass heated by feedback events with
common supernova feedback scheme adopted by all mesh-based
other particle-based codes. Gas particles that receive feedback
codes, we refer the readers to Section 3.2.
compute their required CFL timestep at the timestep prior to
receiving feedback, which helps to prevent numerical instabil-
5.5. CHANGA ity and overcooling. GRACKLE cooling is implemented by
CHANGA is a reimplimentation of GASOLINE (see Section 5.6) applying a half timestep of cooling, then a full timestep of
in the CHARM++ runtime system.61 CHARM++ (Kale & external PdV heating, followed by a final half timestep of
Krishnan1996,pp.175–213)62enablestheoverlapofcomputation cooling, as in CHANGA (Section 5.5).
and communication and provides adaptive load balancing
infrastructure,allowingCHANGAtoscaletohundredsofthousands
of processor cores (Menon et al. 2015). The hydrodynamics in 5.7. GADGET-3
CHANGA closely follows that of GASOLINE. SPH forces are GADGET-3 is an updated version of GADGET-2, a cosmo-
calculated using the method of Ritchie & Thomas (2001), and logical tree-particle-mesh SPH code that was originally
energyisdiffusedusingtheschemeofShenetal.(2010),bothof developed by V. Springel (Springel et al. 2001; Springel
which providing a more accurate treatment of multi-phase ISM. 2005).65 GADGET-3 has important updates from GADGET-2,
Timestepsaredeterminedbytheminimumofanaccelerationand
suchasdomaindecompositionanddynamictreereconstruction
a CFLcriterion.Furthermore,the timesteps of neighbors arekept
whichmayslightlyaltertheN-bodydynamics.TheGADGET-3
withinafactorof2ofeachotherasinSaitoh&Makino(2009)in code used in this comparison is a modified version of the
order to accurately integrate highly supersonicflows.
original GADGET-3 by K. Nagamine and his collaborators,
For this work, a kth nearest-neighbor algorithm is used to
which includes pressure-entropy formulation by Hopkins
find the Nsmooth = 64 nearest neighbors which are smoothed (2013), time-dependent artificial viscosity, variable smoothing
with the Wendland C4 kernel (Dehnen & Aly 2012) to lengths, among others (e.g., Choi & Nagamine 2012; Thomp-
determine hydrodynamic properties. Unlike conventional son et al. 2014; Aoyama et al. 2016).
versions of CHANGA or GASOLINE, the supernovae thermal
For thepresent study,the GADGET-3group adopts a quintic
energy, mass, and metals are directly distributed to the 64 spline smoothing kernel (Morris 1996) with N = 64. The
neighboring gas particles.63 Gas particles that are neighbors of ngb
implementation of supernova feedback is based on an updated
particlesthatwillexplodeasasupernovaintheirnexttimestep
versionofTodoroki(2014)thatlargelyfollowsaSedov–Taylor
are put on timesteps suitable for their post-supernova thermal
blast wave method outlined in Stinson et al. (2006, 2013, but
energy, preventing them from being on a much smaller
not their cooling shutoff model). The exact model used in the
timestep required in the CFL condition. GRACKLE cooling is current work isfully described inAoyama etal.(2016)but, in
implemented but it does not self-consistently account for the
brief,theimplementationcomprisesthefollowingsteps.Every
PdV work or other external sources of energy, a requirement
for CHANGA and GASOLINEʼs energy integration. Therefore, time a star particle explodes, we compute the “shock radius”
based on Chevalier (1974) and McKee & Ostriker (1977), and
wesplittheenergyintegrationintoahalftimestepofGRACKLE
cooling, then a full timestep of PdV heating, and finally a then find the gas particles within the radius. We then inject
thermalenergyandmetalyieldsintotheidentifiedgasparticles
second half timestep of cooling.
within the shock radius, weighted by the SPH spline kernel.
60 Thismeansthatinordertoretrievegasinternalenergyortemperature(e.g., Finally, we note that the results of this version of GADGET-3
Section6),thepressuresupporttermneedstobesubtractedoutfromRAMSESʼs are not representative of all the GADGET-3 codes in the
pressurefield,whichistheonlyfieldbeingtracked. community,becausesomeoftheresultsarestronglydependent
61 Thewebsiteishttp://www-hpcc.astro.washington.edu/tools/changa.html. on the detailed implementations of baryonic physics, such as
62 Thewebsiteishttp://charm.cs.uiuc.edu/. star formation and feedback.
63 Thefeedbackprescriptionusedinthisexperimentneededareimplementa-
Stiotinnsoofntheet faele.d2b0a0c6k).roIuntinpearntiocrumlaarl,lyinusperdeviniouCsHAwNoGrkA,atnhde GsuApSeOrnLoINvEae(er.agt.e, 64 NotethedifferenceinNsmoothfromCHANGAinSection5.5.Dehnen&Aly
determinedbythestellarageandIMFisconvertedtoanenergyinjection“rate” (2012) showed that this kernel can use larger neighbor numbers without the
whichisthenincorporatedintothethermalenergyintegrationofneighboring pairinginstabilitywhichmayeffectivelyremoveresolution,andthatdoingthis
gas particles. By contrast, for this study, a supernova event occurs improvesperformanceonanumberofbasichydrodynamicstests.
instantaneously,makingtherateanill-definedquantityintheirexistingenergy 65 The website is http://www.h-its.org/tap-software-en/gadget-code/ or
integrationmachinery. http://www.mpa-garching.mpg.de/gadget/.
7
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
Figure2.500MyrcompositeofgassurfacedensitiesfromSim-noSFwithradiativegascoolingbutwithoutstarformationorsupernovafeedback.Eachframeis
centeredonthegalacticcenter—locationofmaximumgasdensitywithin1kpcfromthecenterofgasmass.Forvisualizationsoftheparticle-basedcodeshereafter
(Figures1–3,14–15,32,34,and35)—butnotinanyotheranalysesexceptthesefigures—ytusesanin-memoryoctreeonwhichgasparticlesaredepositedusing
smoothingkernels.SeeSection5fordescriptionsofparticipatingcodesinthiscomparison,andSection6.1foradetailedexplanationofthisfigure.Comparewith
Figure 14. Simulations performed by: Daniel Ceverino (ART-I), Robert Feldmann (ART-II), Mike Butler (ENZO), Romain Teyssier (RAMSES), Spencer Wallace
(CHANGA),BenKeller(GASOLINE),Jun-HwanChoi(GADGET-3),YvesRevaz(GEAR),andAlessandroLupi(GIZMO).Thehigh-resolutionversionsofthisfigureand
articleareavailableattheProjectwebsite,http://www.AGORAsimulations.org/.
Figure3.SameasFigure2butforSim-SFFwithstarformationandfeedback.SeeSection6.1foradetailedexplanationofthisfigure.SeealsoSection3.2forthe
commonstarformationprescriptionandtheguidelineforsupernovafeedback,andSection5fortheexactdepositschemeofthermalfeedbackenergyimplementedin
eachcode.ComparewithFigures15,21,29,32,34,and35.
5.8. GEAR 5.9. GIZMO
GEAR is a self-consistent, fully parallelized, chemo-dynamical GIZMO (Hopkins 2015) is a new mesh-free Godunov code
tree SPH code (Revaz & Jablonka 2012) which is built on the based on discrete tracers, aimed at capturing the advantages of
publicly available GADGET-2 code (see Section 5.7; Springel bothLagrangianandEuleriantechniques.66Thenumericalscheme
2005). The simulations reported here are run with the improve- implemented in GIZMO, initially proposed by Lanson & Vila
ments dicussed in Revaz et al. (2016), including the pressure- (2008),followstheimplementationofGaburov&Nitadori(2011)
entropyformulationproposedbyHopkins (2013),individualand and relies on the discretization of the Euler equations of
adaptive time-stepping schemes (Durier & Dalla Vecchia 2012), hydrodynamics among a set of discrete tracers. Unlike in the
artificial viscosity (Monaghan & Gingold 1983) supplemented moving mesh technique, where the volume is partitioned by a
withtheBalsaraswitch(f fromBalsara1995),andparticle-based Voronoi tessellation, GIZMO distributes the volume fraction
ij
time-dependent viscosity coefficient (Rosswog et al. 2000). assigned to the tracers through a kernel function. For the current
For this study, the standard cubic spline smoothing kernel work, GADGETʼs standard cubic spline smoothing kernel is used
(Monaghan & Lattanzio 1985) in GADGET-2 is used with with Nngb = 32. Note, unlike SPH codes, these tracers only
N = 50.Thefeedbackenergy,mass,andmetalinjectioninto represent unstructured cells, sharing an “effective face” with the
thnegbISM is implemented following the standard SPH scheme. neighboring cells.67 The Riemann problem is then solved across
these faces using a Godunov method as in mesh-based codes, to
Theimplementationcomprisesthefollowingsteps.Everytime
a star particle explodes, we first find the nearest gas particles, accurately resolve shocks without artificial dissipation terms.
according to the weighted number of neighbors as defined in Unlike mesh-based codes, these cells are not fixed in space and
Springel & Hernquist (2002). A desired number of neighbors time, resulting in the scheme’s Lagrangian behavior with
N = 50 is used. Then we inject thermal energy and yields
ngb 66 Thewebsiteishttp://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html.
into the neighboringgas particles, weighted by theSPHspline
kernel.GRACKLEcoolingisperformedafterthekickstep,once “67smItooitshinwgorltehngntho”tinign SthPaHt)thdeoeksenrnoetlpslaizyeaninyGroIlZeMiOn t(hwehdaytniasmcicasl,lebdutthies
gas particles have eventually received supernova feedback simply related to each cell’s “effective volume.” For N =32, a radius
ngb
energy and once the size of the next timestep is known. The enclosingapproximatelythismanyneighborsisusedtoestimatetheeffective
adiabaticcooling/heatingisfirstapplied,andthentheradiative vvoolluummee”pceormpinagrtifcrloemattheserceogniodn-owrditehrinacacsuirnagclye,inwteitrh-pamrtoicslteosfeptahreati“oenffleecntgivthe
one provided by GRACKLE. aroundthetracer.
8
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
Figure4.Cylindricallybinnedgassurfacedensityprofilesat500MyrforSim-noSFwithoutstarformationorfeedback.Thecylindricalradiusisfromthegalacticcenter
(locationofmaximumgasdensitywithin1kpcfromthecenterofgasmass;inallanalysesforparticle-basedcodeshereafter)exceptthegraphicalvisualizationssuchas
Figures1–3—rawparticlefieldsareused,nottheinterpolatedorsmoothedfieldsconstructedinyt.Showninthebottompanelisthefractionaldeviationfromthemean
oftheseprofiles.SeeSection6.1formoreinformationonthisfigure.They-axisrangeofthetoppaneliskeptidenticalamongFigures4–7and22foreasiercomparison.
intrinsically adaptive resolution. When the time evolution of the 500Myr snapshots of the participating simulations in two
commonfacebetweentwocellsisconsidered,weusethesecond- setups: Sim-noSF and Sim-SFF. As defined in Section 3, Sim-
order accurate Meshless Finite Mass method (described in noSFreferstoarunwithradiativegascoolingbutwithoutstar
Hopkins 2015). GIZMOʼs time-stepping scheme is fully adaptive, formation or feedback, and Sim-SFF refers to a run with
andcloselyfollowsGADGET-3orAREPO(Springel2010).Italso radiative cooling, star formation and supernova feedback.
includes a timestep limiter by Saitoh & Makino (2009) (see In our simulation analyses, a key role has been played by the
Section 5.5). AGORA recommended community-driven analysis platform yt
GIZMOʼs gravity solver is based on the tree algorithm (Turk et al. 2011; Turk & Smith 2011; Turk 2013, see footnote
inherited from GADGET-3, itself descending from GADGET-2 43). It natively processes data from all nine participating
(see Section 5.7; Springel 2005). Gravitational softenings in simulationcodesdiscussedinthispaper,plusmanyothermodern
GIZMO can be fixed or fully adaptive, but in the reported runs astrophysics codes such as ATHENA (Stone et al. 2008), FLASH
fixed softening length is used matching SPH codes. To model (Fryxell et al. 2000), GADGET-3-SPHS (Read & Hayfield 2012),
supernovafeedback,theGIZMOsimulationsshownhereadopts NYX(Almgrenetal.2013),andORION(Trueloveetal.1998),to
asimilarfeedbackstrategyusedinGEAR(Section5.8).Thatis, name a few. Interested readers may try a unified, publicly
energy,massandmetalsaredistributedamongtheneighboring availableytscriptemployedinthepresentanalysesthathasbeen
gas particles/cells in a kernel-weighted fashion, but with developedthroughouttheprogressofthisstudy.68,69Wealsoplan
N = 32. For star particles, timesteps are constrained to to make data sets used in the present study publicly available in
ngb
prevent supernovae from exploding in the timestep when the the near future (see Section 7 for more information).
stars formed. Lastly, we caution that different feedback
implementations using GIZMO in the literature adopt different 6.1. Gas Disk Morphology
algorithms to distribute supernova feedback energy (e.g.,
Hopkins et al. 2014, by the FIRE Collaboration). Wefirstexaminethemorphologyofgasdisksevolvedineach
ofthecodesinourexperiments.InFigures1–3,wecompilenine
panels that exhibit the results of the isolated disk galaxy
6. RESULTS
Inthissection,welayouttheresultsofthefirstisolateddisk 68 Thewebsiteishttp://bitbucket.org/mornkr/agora-analysis-script/.
69 ytversion3.3orlaterisrequiredforthescripttoreproduceouranalyses.
galaxycomparisonbytheAGORACollaboration.Wefocuson ForthefiguresandplotsinSection6,theyt-devchangesetd7f213e1752e
similarities and discrepancies discovered by comparing isused.
9
TheAstrophysicalJournal,833:202(34pp),2016December20 Kimetal.
Figure5.SameasFigure4butforSim-SFFwithstarformationandfeedback.ComparewithFigures22and27.
simulations, first with radiative gas cooling but without star be seen in this figure. To model the gas disk, SPH particles or
formation or feedback, Sim-noSF, and second with star GIZMOʼs discrete tracers are generated by drawing random
formation and feedback, Sim-SFF, by the nine participating numbers from the distribution function given by an analytic
codes. Each panel displays the disk gas surface density in a densityprofile.ThisbydefinitionresultsinPoissonnoiseinthe
30kpc box centered on the location of maximum gas density disk surface density shown here. Readers can also observe
within1kpcfromthecenterofgasmass.Thiscenteringcriterion slightdifferencesbetweenmesh-basedandparticle-basedcodes
isadoptedinallsubsequentfiguresandplots.Forvisualizations inhowthedensity fieldisrepresented intheircalculations.By
oftheparticle-basedcodeshereafter(Figures1–3,14–15,32,34 the nature of reconstructing the density from the positions of
and35)—butnotinanyotheranalysesexceptthesefigures—yt particles, the particle-based codes may smooth out the strong
uses an in-memory octree on which gas particles are deposited density contrast in the IC at the edge of the initial gas disk.
using smoothing kernels. The resolution of this octree governs Interestingly, in Figures 2 and 3 at 500Myr, there are other
the resolution of produced images. If more than eight particles subtledifferencesnoticeablebetweenmesh-basedandparticle-
are in an oct, that oct is refined into 64 child octs (i.e., yt basedcodesaswellaswithinthesesub-groups.Whilethepeak
parameters n_ref = 8 and over_refine_factor = 2), densitiesandfilamentarystructuresinthedisksareverysimilar
providing compatible or better image resolution than a typical across all codes, it is noticeable that the mesh-based codes
SPH visualization. The densities are assigned to the octree in a typically show lower densities in the inter-arm regions of the
scatter step. That is, we first calculate a particle’s smoothing disks. The typical densities in those inter-arm regions—while
length,andthenaddtheparticle’sdensitycontributiontoallcell notcontainingmuchmass—maydifferbyasmuchasanorder
centersoftheoctreecellsthatarewithintheparticle’ssmoothing of magnitude between mesh-based and particle-based codes
sphere. (see also Section 6.7 and Figure 35 for a related discussion on
We asked every code to output the state of the simulation spatial resolution). Another distinguishing aspect among the
immediately after it was initialized—so-called “0Myr snap- participatingcodesisthenumberofdenseclumpsformed.This
shot”—to allow ourselves to directly compare whether the IC is true in both the simulations with and without star formation
generation was successful and consistent among codes. This andfeedback(seealsoSection6.4andFigures21and23fora
exercise has been strenuously carried out for all the analyses related comparison of newly formed stellar clumps).
items presented in Sections 6.1–6.3 and 6.7, enabling us to Figures 4 and 5 are the cylindrically binned gas surface
correct inconsistently initialized simulations early inthe study. density profiles for Sim-noSF and Sim-SFF, respectively. The
One such example, the surface density comparison of 0Myr cylindrical radius is defined as the distance from the galactic
snapshots,isshowninFigure1.Acleardistinctioningasdisk center. Raw particle fields are used for profiles of the particle-
initializationbetweenmesh-basedandparticle-basedcodescan based codes, not the interpolated or smoothed fields
10