A new paradigm for reproducing and analyzing N-body simulations of planetary systems

0.93 MB·

by  Hanno Rein

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed26January2017 (MNLATEXstylefilev2.2)
A new paradigm for reproducing and analyzing N-body simulations
of planetary systems
Hanno Rein1,2 and Daniel Tamayo1,3,4
1DepartmentofPhysicalandEnvironmentalSciences,UniversityofTorontoatScarborough,Toronto,OntarioM1C1A4,Canada
2DepartmentofAstronomyandAstrophysics,UniversityofToronto,Toronto,Ontario,M5S3H4,Canada
7 3CanadianInstituteforTheoreticalAstrophysics,60St.GeorgeSt,UniversityofToronto,Toronto,OntarioM5S3H8,Canada
1 4CentreforPlanetarySciencesFellow
0
2

n Accepted2017January24.Received2017January20;inoriginalform2016November18
a
J

5 ABSTRACT
2 The reproducibility of experiments is one of the main principles of the scientific method.
However,numericalN-bodyexperiments,especiallythoseofplanetarysystems,arecurrently

] not reproducible. In the most optimistic scenario, they can only be replicated in an approx-
P
imate or statistical sense. Even if authors share their full source code and initial conditions,
E differences in compilers, libraries, operating systems or hardware often lead to qualitatively
h. differentresults.
p We provide a new set of easy-to-use, open-source tools that address the above issues,
- allowingforexact(bit-by-bit)reproducibilityofN-bodyexperiments.Inadditiontogenerat-
o ingcompletelyreproducibleintegrations,weshowthatourframeworkalsooffersnoveland
r
t innovativewaystoanalyzethesesimulations.Asanexample,wepresentahigh-accuracyin-
s
tegrationoftheSolarSystemspanning10Gyrs,requiringseveralweekstorunonamodern
a
[ CPU. In our framework we can not only easily access simulation data at predefined inter-

valsforwhichwesavesnapshots,butatanytimeduringtheintegration.Weachievethisby
1
integrating an on-demand reconstructed simulation forward in time from the nearest snap-
v
shot.Thisallowsustoextractarbitraryquantitiesatanypointinthesavedsimulationexactly
3
(bit-by-bit),andwithinsecondsratherthanweeks.
2
4 WebelievethatthetoolswepresentinthispaperofferanewparadigmforhowN-body
7 simulationsarerun,analyzed,andsharedacrossthecommunity.
0
Keywords: methods:numerical—gravitation—planetsandsatellites:dynamicalevolution
.
1 andstability
0
7
1
:
v
1 INTRODUCTION arenotindividuallyreproducible.Thecentralproblemisthatmost
i
X dynamicalsystemsofinterest(e.g.,theSolarSystem)arechaotic.
Scientists have studied the motion of gravitationally interacting
r Thus,anydiscrepancyintherepresentationoffloatingpointnum-
a bodies such as the planets in the Solar System for centuries (for bers amplifies exponentially fast. This means that seemingly in-
a comprehensive historical review see Laskar 2012). Today, such significantdifferencesinhardware,software,initialconditions,or
studies are mostly done on computers through N-body experi-
even the choice of which times to output data can change, e.g.,
ments. Unsurprisingly, answering new questions often requires
whethertwoplanetscollideornotinagivensimulation.
pushing the limits of current technology. When studying the dy-
namicalevolutionofastrophysicalsystems,thisoftentranslatesto For example, Laskar & Gastineau (2009) showed by direct
severalmonthsofcomputingtime. integration that in about 1% of their simulations, Mercury col-
Working with such simulations can be cumbersome. Tradi- lideswithanotherterrestrialplanetorplungesintotheSuninthe
tionally, one has to specify both the quantities to record and the next5Gyr.Sincethen,severalothergroupshavestudiedthephys-
output cadence ahead of time. But this is not easy to do, as the icaloriginofthisinstability(Lithwick&Wu2011;Batyginetal.
appropriatechoicesareoftenonlyclearafteraninitialanalysisof 2015),butcouldnotdirectlycomparetheirtheorytothe6.2million
the data. The problem is compounded if one considers that other CPU-hourdatasetproducedbyLaskar&Gastineau(2009).Aswe
groups (perhaps with different goals) might also want to analyze explainfurtherbelow,ifoneisinterestedinanalyzingthesimula-
theoutputdataatalatertime. tionsresultingincollision,itisnotenoughtohavethesamecode
Evenworse,currentN-bodyintegratorsyieldtrajectoriesthat andexactinitialconditions.Onewouldhavetoexpendcomparable
(cid:13)c 0000RAS
2 HannoRein,DanielTamayo
computational time to generate many simulations, and select the access to the actual source code used in a simulation, there is no
(different!)onesthatleadtocollisions. hopeofachievingbit-wisereproducibility.
Chaos in N-body integrations therefore not only leads to in- Unfortunately,itisoftenoverlookedthateventheexactsource
efficient use of computational resources, it is also a fundamental code(orevenbinaryexecutable)andinitialconditionsmightpro-
challengetoreproducibilityandabarriertothescientificprocess. ducedistinctoutcomesondifferentmachines.Machine-dependent
Inthispaperwepresentsolutionstothetechnicalobstaclesham- resultscanoriginatefromphysicaldifferencesinhardware,compil-
pering reproducibility, and propose a new paradigm for perform- ers,compileroptions,operatingsystemsorlibraries.Infact,onlya
ingN-bodyintegrations.Inourframework,machine-independent verylimitedsetoffloatingpointoperationssuchasadditionsand
integrationalgorithmsgeneratereloadablebinarysnapshotsofthe multiplicationsareguaranteedbytheIEEE754floatingpointstan-
simulation as it progresses. The resulting Simulation Archive can dard(seeISO2011)togivethesameresultsonallmachinesthat
then be shared across the community, and different members can comply with the standard. We note in particular that often-used
robustlychooseafterthefactwhatdynamicalquantities(e.g.,or- functions such as sin(), cos(), and pow(), and thus any codes
bitalelements)toextract.Animportantdifferenceinthisdataanal- thatemploysthem,arenotmachineindependent.
ysisstageisthat(unlikewiththeoriginalsimulation)thetaskcan Finally,oneobviouslyrequiresasetofinitialconditionstore-
beeasilyparallelized.Wehaveimplementedthesefeaturesinthe produce an N-body simulation. Bit-wise reproducibility addition-
open-sourceN-bodypackageREBOUND. allydemandsstoringtheseinitialconditionsinabinaryformat,in
Although, for the sake of definiteness, we focus on the So- ordertohavetheexactfloatingpointrepresentationofallnumbers;
larSystemandtheHLTausysteminthispaper,ourdiscussionis atext-basedformatisnotsufficient.
equallyapplicabletoawidevarietyofdynamicalsystems.Thepa-
perisorganizedasfollows.Wefirstdiscussindetailtheissueof
reproducibilityinSec.2.Then,wedescribetheconceptoftheSim-
3 THESIMULATIONARCHIVE
ulationArchiveandhowitcanbeusedtoefficientlyanalyzelong-
termsimulationsinSec.3.WelookattwoexamplesinSec.4,a Wenowdescribehowtoovercomethebarrierstoreproducibility
longtermintegrationoftheSolarSystem(Sec.4.1)andanunsta- mentioned above. First, it is essential to work with open-source
blesystemresemblingtheorbitalparametersoftheHLTausystem software and tools. We work with REBOUND (Rein & Liu 2012),
(Sec.4.2).WeconcludewithasummaryinSec.5. and in particular the integrators WHFast (Rein & Tamayo 2015)
andIAS15(Rein&Spiegel2015).WHFastisaheavilyoptimized
implementationofthestandardsecond-orderWisdom-Holmanal-
2 BARRIERSTOREPRODUCIBILITY gorithm(Kinoshitaetal.1991;Wisdom&Holman1991)andcan
beusedwithsymplecticcorrectorsuptoorder11(Wisdom2006)
Current N-body simulations, especially those of chaotic systems, inbothJacobiandheliocentriccoordinates.IAS15isahigh-order,
are only reproducible in a statistical sense. Different groups em- adaptive Gauss-Radau scheme that is accurate to machine preci-
ployingthesamesimulationwillonlyagreeonthedistributionsof sionandcanhandlebothconservativeandnon-conservativeforces.
outcomes from a large enough number of simulations. However, BoththesymplecticWHFastandnon-symplecticIAS15integrators
oneisofteninterestedinparticularoutcomes,e.g.,Mercurydiffus- aremachineindependent.
ing onto a collisional trajectory with Venus (Laskar & Gastineau Weachievethisbysubstitutingfunctioncallstomathematical
2009), or Jupiter ejecting a fifth giant planet from the Solar Sys- librarieswithappropriateseriesexpansions.Forexample,weim-
temtoexplainthecurrentorbitalarchitectureofthegiantplanets plementaseriesexpansionofStumpfffunctionsforWHFast(Rein
(Nesvorny`&Morbidelli2012).Insuchacase,twogroupswillnot & Tamayo 2015). For IAS15, we implement our own routine to
beabletoreproducethesameoutcomeeveniftheysharetheirex- calculate7throotsinordertoavoidcallingthemachine-dependent
act initial conditions. To obtain these rare cases, one would have pow()function.Implementingourownmathematicalfunctionsof-
toexpendthecomputationalcostofalargesuiteofsimulationsto ten also leads to a more accurate and faster algorithm (Rein &
extractafewqualitativelysimilarbutdistincttrajectories.Thisis Tamayo2015).WealsonotethatREBOUNDiswritteninC99,apro-
notonlygrosslyinefficient,butalsorendersdetailedcomparisons grammingstandardthatexplicitlystatesthatfloating-pointopera-
impossible. tionsmaynotberearrangedbythecompilerunlessitcanguarantee
Oneofthefirstobstaclesinreproducibilityofnumericalsimu- thattheresultisbit-wiseequivalent.Thisisstilltrueifoptimization
lationsisthelackofaccesstosourcecode(Baker2016).Although flagssuchas-O3areused,butnotif--fast-mathisenabled.
manyauthorshavepublishedtheirnumericalmethodsnotonlyas Thus,byusingREBOUNDandoneoftheaboveintegrators,sim-
papers but with the accompanying code (Chambers & Migliorini ulationsareautomaticallymachineindependent.Ifauthorschoose
1997;Duncanetal.1998;Rein&Tamayo2015),othershavenot. to share their initial conditions, simulations are then in principle
Whereas one should in principle be able to reconstruct the algo- bit-wise reproducible. However, several additional technical de-
rithmfromadescriptioninapaper,thisisoftennotpractical,and tails haveto be overcomein orderto replicate thesame series of
iseffectivelyimpossiblewhenthegoalisbit-wisereproducibility. timesteps.Inthefollowingsectionwedescribeourconceptofthe
Asasimpleexample,considertheexpressions SimulationArchive,whichaddressestheseissues,andallowsusers
y =(2/3)·x and y =(2·x)/3 torecreateandrestartintegrationsatarbitrarytimesalmostinstan-
1 2
taneously. This enables exciting new ways to analyze and share
whicharemathematicallyequivalentbutwillyielddifferentresults simulationdata.
onacomputerwithfinitefloatingpointprecision1.Thus,without
Bycontrast,y1 mightleadtobiasedresultsbecause2/3evaluatestothe
1 Notethaty2isingeneralthebetterchoicebecauseitrandomizeserrors. sameroundednumbereachtimetheexpressionisevaluated.
(cid:13)c 0000RAS,MNRAS000,000–000
AnewparadigmforreproducingandanalyzingN-bodysimulationsofplanetarysystems 3
3.1 Asinglebinaryfile of mixed-variable symplectic integrators (e.g., Saha & Tremaine
1992),whicharethedefactostandardforlong-termplanetaryin-
General purpose binary file formats for N-body simulations have
tegrations,tonon-symplecticintegrators,andtotheincorporation
alreadybeendevelopedbyseveralothergroups(seee.g.Faberetal. ofadditionaleffectsbeyondpoint-sourcegravitationalforces.
2010; Farr et al. 2012; Cai et al. 2015). This paper is not about
anewbinaryformatfor N-bodysimulations.Rather,wedescribe
howtostoredatainawaythatallowssimulationstoberestarted
in a way that reproduces the original simulation exactly. To our 3.3 Mixed-variablesymplecticintegratorsand
knowledge,thishasneverbeendonebeforeinapubliclyavailable synchronization
codethatisusedforlongtermplanetarysimulations.
When trying to restart a simulation using a mixed-variable sym-
WhatwerefertoastheSimulationArchiveisasinglebinary
plectic (MVS) integrator like WHFast, one encounters three syn-
filethatincludesallrelevantinformationtoreproducean N-body
chronizationissues.Wedetailhowweovercomethesechallenges
simulationdowntothelastbit.Inparticular,westore:
inWHFastspecifically,butemphasizethatthesepitfallsandsolu-
- Nameandversionoftheprogramusedtogeneratethefile tionsareapplicabletoMVSintegratorsingeneral.
- Initialconditions First, MVS integrators split the evolution of the system into
- Allconstantsandconfigurationparameters,suchasgravitational exactly solvable steps, which are then interleaved with one an-
constantandtimestep other.TheWHFastintegrator,likemanyotherMVSintegrators,is
aleapfrog-styledrift-kick-driftscheme.Thedriftstepcorresponds
Then,asthesimulationprogresses,weappendthebinaryfilewith toanoperatorevolvingthesystemundertheKeplerianpartofthe
snapshotsofalldataneededtorestartthesimulationatthattime. Hamiltonianandthekickstepevolvesthesystemundertheinter-
In the case of a simulation of the Solar System with eight plan- actionterm.Eachofthetwodriftstepsevolvesthesystemforhalf
ets,wetypicallyappendasnapshotevery50000yrstotheSimu- atimestep.Inalongsimulation,onecanspeedupthecalculation
lationArchive.Eachsnapshotisonly500byteslongbutcontains bycombiningthelastandfirstdriftstepsofadjacenttimestepsinto
allthenecessaryinformationtoexactlyrestartasimulationatthat one.However,ifpositionsandvelocitiesarerequiredafterapar-
time.Fora10Gyrintegration,thetotalfilesizeoftheSimulation ticulartimestep,onemustreintroducethehalf-timestepdriftstep,
Archiveisabout100MB.Cosmologicalorotherlarge-Nsimula- aprocessthatwerefertoassynchronizingthepositionsandveloc-
tionswouldhavetofindadifferentbalancebetweenfilesizeand ities.Thisleadstoaproblemifoneconsiderstwointegrationsof
thecomputationtimebetweensnapshots. thesametrajectoryacrossaparticulartimestep,whereinonecase
We note that our implementation does currently not support onechoosestooutputparticleinformation(andthusmustsynchro-
big-endianmachines.However,sinceallmodernx86/x64architec- nize),andintheotheronedoesnot(andthusperformsafulldrift
turesarelittleendian,weseefewreasonstoimplementthisfeature. step). While evolving the system under the drift operator for two
Shouldtheneedarise,itisstraightforwardtoextendtheformatto half-timesteps(inordertooutputdata)ismathematicallyequiva-
beendianagnostic. lenttoafulldriftstep,thetwowillgiveadifferentnumericalresult
Furthermore,wenotethatourimplementationsupportsboth becauseofthefinitefloatingpointprecision.Todealwiththis,the
32 and 64 bit architectures. However, files created on 32 bit ma- SimulationArchivealwaysstorestheparticlepositionsandveloc-
chinescancurrentlynotbereadona64bitmachineandviceversa. itiesintheunsynchronizedstatefollowingthekickstep.Ifanad-
Giventhewidespreadadoptionof64bitsystemsthesedays,wesee ditionalhalfdriftstepisrequiredtosynchronizethepositionsand
fewreasonstoimplementtheSimulationArchiveinanarchitecture velocities, the result is then discarded before the integration con-
agnosticway.Shouldtheneedarise,thisfeatureisstraightforward tinueswithafulldriftstepfromthestatethatwascachedbythe
toimplementanddoesnotaffectthereproducibilityofsimulations. SimulationArchive.
Second,thedriftandkickstepsareoftensolvedindifferent
coordinatesystems(e.g.,Wisdom&Holman1991).TheKeplerian
3.2 Extractingarbitraryquantitiesatanytime partisbestsolvedinJacobicoordinates2,whereastheinteraction
Hamiltonianissolvedinaninertialframe.Wethushavetoconvert
ManyN-bodycodesallowonetosavebinarysnapshotsofsimu-
backandforthbetweentwodifferentcoordinatesystems,whichcan
lations.TheSimulationArchivegoesbeyondthisbyguaranteeing
thatdifferentmachinesreloadingandintegratingthesamesnapshot harm precision in long term integrations (Rein & Tamayo 2015).
However, at closer inspection one finds that the interaction step
willgiveidenticalresultsbit-by-bit.
onlyneedstoaccessthepositionsintheinertialframe,nottheve-
Thestoragerequirementsforonesnapshotareverysmall,only
locities. Furthermore, only the accelerations calculated in the in-
500bytesforaWHFastsimulationoftheSolarSystem.Thisallows
ertial frame are needed in Jacobi coordinates. Thus, a simulation
ustostoremanysnapshotswhicharetypicallyspacedonlyafew
can stay in Jacobi coordinates for the entire integration and only
secondsapartinwalltime.Onecanthereforeaccessanyvalueof
calculatethepositionsintheinertialframefortheinteractionstep
any particle at any timestep by loading the snapshot just prior to
(butnotrecalculatetheJacobicoordinatesfromthosepositionsas
therequestedtimeandthenintegratingthesimulationforwardin
positionsdonotchangeintheinteractionstep).IntheSimulation
time.Sincerestoringasimulationfromabinaryfileisalmostin-
Archive,westoretheJacobicoordinates,notthecoordinatesinthe
stantaneousandsnapshotsareonlyafewsecondsapart,onehasto
waitatmostafewsecondsforthisshortintegrationtofinishand
yieldthedesiredquantities.Furthermore,wenotethatthisprocess
iseasilyparallelizableifonewantstoaccessparametersatmultiple 2 Jacobi coordinates are advantageous for systems with planets on well
times(SeeSec.4.1and4.2). separatedorbitssuchastheSolarSystem.Forothersystemswhereorbits
We now address technical details pertaining to the class mightcross,heliocentriccoordinatesmightbebeneficial.
(cid:13)c 0000RAS,MNRAS000,000–000
4 HannoRein,DanielTamayo
inertialframe.Werewetostoretheinertialcoordinatesandthenre- perform a synchronization step. This can lead to a further issue
calculatetheJacobicoordinateslater,wewould,onceagain,geta ifnon-gravitational effectsare present,whichisoften thecasein
differentnumericalresultbecauseoffinitefloatingpointprecision. N-body simulations, for example to model the effects of general
Thus,wewouldnotbeabletorestartsimulationsexactly. relativityortides.
Third,whensymplecticcorrectorsareused,acorrectorstepis In REBOUND these effects are modelled as either additional
performedatthebeginningoftheintegrationandwheneveranout- forces or so called post-timestep modifications. If the effects are
putisneeded(e.g.,Wisdom2006).Similartothecalculationfrom simplyforces,thentheydonotneedtobeappliedforthesynchro-
Jacobicoordinatestotheinertialframeandback,thiscorrectorstep nizationstepwhichonlyinvolvesadriftstep,butnotakickstep
candegradetheprecisionandspeedinlongtermsimulationsifap- (whereforcesarecalculated).However,ifpost-timestepmodifica-
pliedrepeatedly.Andagainsimilartoabove,wewouldnotbeable tionsareused,thentheseneedtobealsoappliedwhenthesynchro-
torestartasimulationexactlywerewetostorethecoordinatesafter nizationstepsareperformed.Furthermore,ifsymplecticcorrectors
applyingasymplecticcorrector.Thisiswhywhenstoringcoordi- areused,theneitherkindofeffectmustbepresentatthetimeof
nates in the Simulation Archive, we do not apply the symplectic synchronizationbecausesymplecticcorrectorsapplyboththekick
correctorstepsbeforehand.Thecorrectorstepwillbeperformedat anddriftstepsrepeatedly.
alatertime.
Insummary,intheSimulationArchive,westoretheJacobico-
ordinatesofallparticlesinwhatwecallanunsynchronizedstate,
i.e.afterthekickstepandwithoutapplyingsymplecticcorrectors.
Onlythecombinationoftheseconceptsallowsustorestartsimula-
4 EXAMPLES
tionsexactlydowntothelastbit.
As a result, if we want to get any physical meaningful out-
WenowgivetwoconcreteexamplesofhowtheSimulationArchive
put from the Simulation Archive, we need to synchronize the co- canbeusedinpractice,onewiththesymplecticWHFastalgorithm,
ordinates beforehand. To reemphasize, this is done at a time the andtheotherwiththeadaptive,high-accuracyIAS15integrator.
simulationisanalyzed,notwhilethesimulationisrunning.Inthis
synchronizationstep,wefirstperformhalfadriftsteptoadvance
positionsandvelocitiestothesametime.Wethenapplythesym-
plecticcorrectorsandconvertfromJacobitoinertialcoordinates.
Ifwewanttorestartasimulation(anddonotrequireanout- 4.1 SolarSystem
put),thenwedonotneedtoperformthesynchronizationstepsand
As an example, we now integrate the eight major planets of our
simply restart from the Jacobi coordinates in the unsynchronized
SolarSystemforwardintime.Thesamediscussionappliestoother
state.
dynamicalsystemssuchasexoplanetarysystems.
WealsoimplementtheSimulationArchivefortheheliocentric
The shortest timescale in the Solar System that we consider
versionofWHFast.Heliocentriccoordinatescanbebeneficialfor
isthe88-dayorbitalperiodofMercury.Thelongesttimescalewe
system in which orbits are crossing (for a discussion of differnet
areinterestedinisroughlytheageoftheSolarSystem,4.6billion
coordinatesystems,seee.g.Hernandez&Dehnen2016).Thesyn-
years, which also happens to approximate the remaining lifetime
chronization issues described above are exactly the same, except
oftheSunonthemainsequence.Thehugeseparationofscales(11
thatwehaveheliocentriccoordinatesinsteadofJacobicoordinates. ordersofmagnitude)makesthisanextremelydifficultproblemto
Indeed,theaboveconsiderationsneedtobeaddressedtocreateany
calculate,eventhoughthephysicallawsgoverningthemotionof
reproducibleMVSintegrator.
planetshavebeenwellknownsinceNewton(ignoringsmallgen-
eralrelativisticandothercorrections).
3.4 IAS15 Becauseofthis,thedynamicalevolutionoftheSolarSystem
stilloffersmanysurprisesforresearcherstothisday.Infact,only
For integrators that do not rely on coordinate transformations inthelast25yearshavedirectnumericalintegrationsofallSolar
restartingsimulationsbit-by-bitismorestraightforward.However, Systemplanetsbecomefeasible.Thisismainlythankstothedevel-
someissuesmightstillarise,forecampleiftheintegratorsdepends opmentoffastcomputersandnovelintegrationmethods,mostno-
oninformationfromprevioustimestepssuchasapredictedvalues tablymixed-variablesymplecticintegrators(Kinoshitaetal.1991;
forthetimestepandcoordinates.WeimplementedtheSimulation Wisdom&Holman1991).Nevertheless,accuratedirectnumerical
ArchivefortheIAS15integrator(Rein&Spiegel2015).Forthis integrationsoftheeightSolarSystemplanetsoverseveralbillion
specificintegrator,oneneedstostoreseveralinternalarrayswhich years may still require months of computing time per simulation
IAS15usestopredictvaluesforthenextstepandspeedupconver- andseveralmillionCPU-hoursforanensembleofsimulations.At-
gence.BecauseaniterationinIAS15mightconvergetoaslightly temptstosignificantlyspeedupsimulationswiththehelpofparal-
different value depending on the initial guess, we need to store lelizationarefutileastheunderlyingproblemisinherentlysequen-
thesepredictedvaluesateverysnapshottoguaranteebit-wisere- tial(althoughattemptshavebeenmade,seeSahaetal.1997).
producibility.Thisresultsinasinglesnapshotbeingaboutafactor A further complication arises because the Solar System is
of10largerthanforWHFast. chaoticontimescalesofseveralmillionyears,whereasthelife-time
ofthesystemisseveralbillionyears.Onethushastointegratean
ensembleofsimulationstodrawstatisticalconclusionsofthesys-
3.5 Additionalforcesandpost-timestepmodifications
tem’slong-termevolution.ForexampleLaskar&Gastineau(2009)
Asdiscussedabove,whenaccessingaWHFastSimulationArchive requiredmorethan6millionhoursofwalltimeforasetof2500
snapshot to output physically meaningful quantities, we have to simulations,eachspanning5Gyrs.
(cid:13)c 0000RAS,MNRAS000,000–000
AnewparadigmforreproducingandanalyzingN-bodysimulationsofplanetarysystems 5
4.1.1 Initialconditionsandintegrator 10-9
FzoonrsosuyrsitnetmegtroatpiornovoidfethuesSwoiltahriSnyitsitaelmco,nwdeitiuosnes3th.eMNoAreSaAccHuroartie- error10-10
y
ephemeridesareavailable(Fiengaetal.2011),butgiventheother nerg10-11
approximations we make (see below), we feel confident that the e e
v
imnietniatallcpohnydsitiicoanlseffweecutss.eareaccurateenoughtocapturethefunda- relati10-12
We use the symplectic Wisdom-Holman type integrator
10-13
WHFast(Rein&Tamayo2015)mentionedabovewithatimestep 107 108 109 1010
time [yrs]
of6days.Toimprovetheenergyconservation,weuseasymplectic
correctoroforder11(Wisdom2006).Withthesesettings,asimu- Figure1. ThisplotshowstherelativeenergyerrorofourSolarSystem
lationtakesapproximatelyonemonthofwalltimeonanIntelXeon simulationasafunctionoftime.Thebluecurveshowstheslopeexpected
CPU(E5-2697v2,2.70GHz). fromanunbiased t1/2errorgrowth.Ittakes1.9secondstogeneratethedata
forthisplotfromtheSimulationArchiveandthenrendertheplot.
4.1.2 Non-gravitationaleffects
y0.45
ur0.40
Weincludeapost-Newtoniancorrectionforgeneralrelativityinthe erc0.35
M0.30
formofasimplepotential(Nobili&Roxburgh1986), of 0.25
Φ = −3G2mM(cid:12)2 1 (1) ntricity 000...112050
GR c2 r2 ce0.05
ec0.00
foreachplanetwithmassm.IntheaboveequationM isthemass 0 2 4 6 8 10
(cid:12) time [Gyrs]
oftheSun,r istheheliocentricdistanceoftheplanet,andG and
carethegravitationalconstantandthespeedoflight,respectively. Figure2. ThegraylineshowstheeccentricityofMercuryasafunctionof
Thisensuresthatwereproducethecorrectapsidalprecessionfre- time.Becausetheeccentricitybehaviourisdominatedbyshort-termoscil-
lations,wealsoplotamovingaveragewitha1.5-million-yearwindowas
quenciesfortheplanets,inparticularthatofMercury.Weignore
all other non-gravitational effects and the Earth-Moon system is ablackline.Ittakes2.3secondstogeneratethedataforthisplotfromthe
SimulationArchiveandthenrendertheplot.
treatedasoneparticle.
thesnapshotsstoredintheSimulationArchive.Thiskindofanaly-
4.1.3 Simulationarchive
sisisextremelyfast.Readingthebinaryfile,restoringthesimula-
WeusetheSimulationArchiveforourintegration.Individualsnap- tions,synchronizingtheintegrator,andcalculatingtheenergyerror
shotsare51434yearsapart,whichcorrespondstoroughly10sec- foreverysnapshot,followedbyproducingtheplot,takeslessthan
onds in wall time. We integrate the system for 10 billion years, 2seconds.
thusendingupwith194000snapshots.Giventhateachsnapshot In Fig. 2, we plot the eccentricity of Mercury as a function
requires500bytes,thisyieldsafinalSimulationArchivebinaryfile oftime(graycurve).Wealsoover-plottherunningaverageovera
ofroughly100MB. windowof1.5Myrsasablackcurvetofilterouthigh-frequency
Withthissetup,wecannowaccessanarbitrarytimeduringthe oscillations. This plot illustrates that Mercury remains stable in
entire10Gyrintegrationwithin5secondsonaverageandwithinat this particular simulation as expected for ≈ 99% of the simula-
most10seconds. tions (Laskar & Gastineau 2009). Again, we used the Simulation
Archivetoproducethisplotafterthesimulationhadfinished.Read-
ingthebinaryfile,restoringthesimulations,synchronizingthein-
tegrator, calculating the orbital elements, calculating the running
4.1.4 Results
averageandproducingtheplottakesonlyabout2seconds.
WehaverunseveralhundredsimulationsoftheSolarSystemover The first two plots illustrate the usefulness of the Simula-
10Gyrswithslightlydifferentinitialconditions(fractionaldiffer- tion Archive in quickly analyzing and visualizing very long in-
ences of 10−10). We will present the full set of simulations and a tegrations. The next figure illustrates how to extract data at a
discussionofthephysicalresultsinafuturepaper.Inthispaper,we higher cadence than provided by the snapshots in the Simulation
onlyuseonesimulationintheensembletodemonstratetheSimu- Archive.Figure3plotstheeccentricityvariationsoftheEarthover
lationArchive. a 800000 yrs interval, starting at the present day (top panel) and
InFig.1weshowtherelativeenergyerrorofthesystemasa 5billionyearsfromnow(bottompanel).Onecanseelargeoscilla-
functionoftime.Thestraightbluelinecorrespondstoasub-linear tionsonaroughly100000yrtimescale,whichareassociatedwith
errorgrowth∝ t1/2.Wecanseethattheerrorgrowthfollowsthis Milankovitch cycles and changes in the Earth’s climate. Because
t1/2 behaviour for the entire 10 Gyr integration, thus confirming thevariationshappenonatimescalecomparabletotheoutputca-
thattheWHFastintegratorisunbiasedandfollowsBrouwer’slaw denceintheSimulationArchive,wecannotrelyonthesnapshots
(Brouwer1937;Rein&Tamayo2015).Wecreatedthisplotusing intheSimulationArchivetogeneratethisplot.Instead,weloada
snapshotthatisnearbyandthenre-integratethesimulationforward
intime,thistimeoutputtingdatamorefrequently.Theentirepro-
3 http://ssd.jpl.nasa.gov/?horizons cess is automated so that a user can simply request a simulation
(cid:13)c 0000RAS,MNRAS000,000–000
6 HannoRein,DanielTamayo
0.06 planetshadacloseencounterafter∼3.5Myrs.Tofurtheranalyze
h
art0.05 thisencounter,werestartthesimulationjustbeforetheencounter
E
eccentricity of 00000.....0000001234 atcthhnoeidusnbwitineto-rcburwyleda-ibslnleitootrtcehccpeuurrroordeuinnutptcltuiyhbteiblciertaeydpseoteansnrscstieeub.drleeTsswihmtihistuahtliasatthniooeynneol.xytTahcopetroosspausumirbbkelleinccollbyowesacleevadauegisnlee--
0 100000 200000 300000 400000 500000 600000 700000 800000
ableintegrator.
time [yrs]
0.08
Earth00..0067
of 0.05
entricity 000...000234 5WeCprOesNenCtLaUneSwIOpNarSadigmforperformingN-bodysimulationsthat
ecc00..0001 allowsresearcherstoreproduceoneanother’sresultsexactly.This
0 100000 200000 300000 ti4m0e0 0[y0r0s]500000 600000 700000+850e09000 is fundamentally different to what has been done before, where,
inthebestcase,simulationswereonlyreproducibleinastatistical
Figure3. ThetopandbottompanelbothshowtheeccentricityoftheEarth
sense.Webelievethiswillimprovethescientificmethod’sability
overa800000yearperiod.Thetoppanelstartsatthepresentday,while tooperateeffectively.
thebottompanelstarts5billionyearsfromnow.Theoscillationsvisiblein
Weachievebit-by-bitreproducibilitybycarefullyimplement-
theseplotsarerelatedtoMilankovitchcycles.Ittakesonly41secondsto
ing machine-independent integrators and storing all initial condi-
generatethedataforthisplotfromtheSimulationArchiveinparallelwith
24threadsandthenrendertheplot. tionsandparametersinabinaryfile.Becauseourapproachallows
everyone (not just the original authors of a study) to reproduce
particular simulations exactly, this opens up new possibilities for
at a given time, and the Simulation Archive will load the nearby sharingandanalyzingdynamicalsystems.Inparticular,thismakes
snapshotandintegrateittothattime.Notethatthisprocessistriv- it possible for the first time to reanalyse particular trajectories in
iallyparallelizable,allowingustospeeduptheanalysisbymak- chaoticsystemsthatleadtointeresting(andperhapsrare)outcomes
ing use of multi-core desktop machines and clusters. While pro- whichotherwisewouldbeunrecoverable.
ducingFig.3,wemadeuseof24threads,renderingbothplotsin WedescribetheconceptofourSimulationArchive,whichen-
only 41 seconds. The bit-by-bit reproducibility of the Simulation ablestheaboveandprovidesaneasy-to-useinterfaceforextracting
Archivemeansthatwearenolongerrestrictedtoapredefinedout- simulationdataatsavedsnapshots,oratintermediatetimes.This
putcadence;wecanalwaysresamplethesimulationafterthefact makesiteasyandefficienttoanalyzethesimulationretrospectively,
toinvestigateinterestingfeatures. andtooutputarbitraryquantities,exactly(bit-by-bit)asifonehad
storedthemintheoriginalrun.Becauseouralgorithmsaremachine
independent, Simulation Archives can be shared between groups,
4.2 Collisionaltrajectory eveniftheyusedifferentoperatingsystems,compilersorlibraries.
Wehopethattheideaspresentedinthispaperwillbecomea
TheSimulationArchivecanalsobeusedtoanalyzesystemswith
new paradigm for running N-body simulations. This includes the
closeencounters,forexampleifstudyingtheejectionofafifthgiant
expectation that when scientists present N-body simulations, that
planet from the Solar System (Nesvorny` & Morbidelli 2012). As
theywouldsharetheirmachine-independentsourcecode,aswellas
anotherexample,wehereconsideranintegrationoffiveplanetsin
theirexactinitialconditionsandallothernumericalparametersand
the HL Tau system, with semi-major axes at the locations of the
constants.Suchapracticewouldmakeitpossibletoverifyandex-
fivemostprominentgapsinthecircumstellardiskasreportedby
tendoneanother’sresults,renderingthescientificenterprisemore
Brogan et al. (2015). We assign the star one solar mass, and the
planets a mass 3·10−4 times smaller, roughly one Saturn mass, efficientandrobust.
The Simulation Archive has been implemented in the
whichwillquicklyleadtocrossingorbitsanddynamicalinstability
REBOUND code which can be downloaded at https://github.
(Tamayo et al. 2015). The orbits are initialized as circular, with
com/hannorein/rebound. iPython notebooks to reproduce the
aninclinationdrawnfromauniformdistributionbetween0and1
plotsinthispapercanbedownloadedathttps://github.com/
degrees.Theremaininganglesaredrawnrandomlyfromauniform
hannorein/reproducibility-paper.
distribution.
We integrate the system for 10 million years using IAS15,
whichisideallysuitedtostudysuchanunstablesystem,givenits
adaptivetimestep.Forsuchunstablesystems,whereequalintervals
ACKNOWLEDGMENTS
insimulationtimemightcorrespondtovastlydifferentintervalsin
computationtimewhileresolvingparticularlycloseencounters,it ThisresearchhasbeensupportedbytheNSERCDiscoveryGrant
isusefultospecifytheoutputcadenceinwalltime.Thisensures RGPIN-2014-04553.Wethankananonymousreviewerforvalue-
thatintegratingfromonesnapshottothenextalwaystakesanequal ablefeedbackthathelpedusimprovethemanuscript.
amountoftime.Inourcase,wesavesnapshotsat10secondinter-
valsinwalltime,yieldingabinaryof∼3MB.
InFig.4,weplotthesemi-majoraxesofallplanetsassolid
REFERENCES
blacklines,andalsoshowtheradialrangebetweentheirperiastra
andapastrashadedingrey.Theinitialoutputsamplinginthetop Baker,M.2016,Nature
plotiscoarse,butwecaneasilyidentifythatthesecondandthird Batygin,K.,Morbidelli,A.,&Holman,M.J.2015,ApJ,799,120
(cid:13)c 0000RAS,MNRAS000,000–000
AnewparadigmforreproducingandanalyzingN-bodysimulationsofplanetarysystems 7
103
U]
A
axis [102
or
aj
m
mi-
e
s
101
0.0 0.2 0.4 0.6 0.8 1.0
time [yrs] 1e7
60
55
50
U]
xis [A45
a
or 40
aj
m
mi-35
e
s
30
25
20
2800000 3000000 3200000 3400000 3600000 3800000
time [yrs]
Figure4. Thesemi-majoraxes(solidblacklines)ofthefiveplanetsinahypotheticalHLTausystem.Thegrayareascorrespondtotheradialrangebetween
theperiastraandapastraoftheplanets.Planetstwoandthreehaveacloseencounterafter∼3.5Myrs.Toanalysethisspecificcloseencounter,werestartthe
simulationjustbeforethecloseencounterandincreasetheoutputcadencetogeneratethebottomplot.ThankstotheexactreproducibilityoftheSimulation
Archive,thesamecloseencounterisreproduced.
Brogan,C.,Pe´rez,L.,Hunter,T.,Dent,W.,Hales,A.,Hills,R., icsandDynamicalAstronomy,50,59
Corder,S.,Fomalont,E.,Vlahakis,C.,Asaki,Y.,etal.2015,The Laskar,J.2012,ArXive-prints
AstrophysicalJournalLetters,808,L3 Laskar,J.&Gastineau,M.2009,Nat,459,817
Brouwer,D.1937,AJ,46,149 Lithwick,Y.&Wu,Y.2011,ApJ,739,31
Cai,M.X.,Meiron,Y.,Kouwenhoven,M.B.N.,Assmann,P.,& Nesvorny`,D.&Morbidelli,A.2012,TheAstronomicalJournal,
Spurzem,R.2015,ApJS,219,31 144,117
Chambers,J.E.&Migliorini,F.1997,inBulletinoftheAmeri- Nobili,A.&Roxburgh,I.W.1986,inIAUSymposium,Vol.114,
canAstronomicalSociety,Vol.29,AAS/DivisionforPlanetary Relativity in Celestial Mechanics and Astrometry. High Preci-
SciencesMeetingAbstracts#29,1024 sion Dynamical Theories and Observational Verifications, ed.
Duncan,M.J.,Levison,H.F.,&Lee,M.H.1998,AJ,116,2067 J.Kovalevsky&V.A.Brumberg,105–110
Faber,N.T.,Stibbe,D.,PortegiesZwart,S.,McMillan,S.L.W., Rein,H.&Liu,S.-F.2012,A&A,537,A128
&Boily,C.M.2010,MNRAS,401,1898 Rein,H.&Spiegel,D.S.2015,MNRAS,446,1424
Rein,H.&Tamayo,D.2015,MNRAS,452,376
Farr, W. M., Ames, J., Hut, P., Makino, J., McMillan, S., Mu-
Saha,P.,Stadel,J.,&Tremaine,S.1997,AJ,114,409
ranushi, T., Nakamura, K., Nitadori, K., & Portegies Zwart, S.
Saha,P.&Tremaine,S.1992,AJ,104,1633
2012,NewAstron.,17,520
Tamayo, D., Triaud, A. H. M. J., Menou, K., & Rein, H. 2015,
Fienga,A.,Laskar,J.,Kuchynka,P.,Manche,H.,Desvignes,G.,
ApJ,805,100
Gastineau,M.,Cognard,I.,&Theureau,G.2011,CelestialMe-
Wisdom,J.2006,AJ,131,2294
chanicsandDynamicalAstronomy,111,363
Wisdom,J.&Holman,M.1991,AJ,102,1528
Hernandez,D.M.&Dehnen,W.2016,ArXive-prints
ISO. 2011, ISO/IEC/IEEE 60559:2011 Information technology
—MicroprocessorSystems—Floating-Pointarithmetic(New
York,NY,USA:IEEE),58
Kinoshita,H.,Yoshida,H.,&Nakai,H.1991,CelestialMechan-
(cid:13)c 0000RAS,MNRAS000,000–000