I’m at the Numerical Cosmology 2012 workshop in Cambridge this week. Today, the first day, was pretty much dedicated to simulations, focusing mostly on comparisons between different techniques.
Volker Springel was showing off his moving mesh code, AREPO, which uses a Voronoi tessellation scheme to optimally resize mesh cells. The videos on his website are mesmerising. (Neat Voronoi fact: the Voronoi tessellation of a set of points is dual to the Delaunay triangulation.) Why go to all that trouble to make the mesh move around? Well, as was repeatedly pointed out during the day’s talks, Cartesian mesh methods aren’t Galilean invariant (simple translations of objects change their properties, when they shouldn’t), and smooth particle hydrodynamics (SPH) methods suppress (physical) fluid instabilities due to a spurious numerical “surface tension” property that they have. The moving mesh helps to get around these problems – the mesh follows the particle velocities, which ends up guaranteeing Galilean invariance. The use of Voronoi diagrams to define the mesh turns out not to add too much computational expense (efficient algorithms exist), but it does have problems with uniqueness of tessellations (e.g. there’s more than one valid way to tessellate a 4-point square), and solving the Riemann problem for irregularly-shaped cells seems to be harder too. Still, it looks neat.
Volker has apparently been quite critical of SPH methods in recent years – the method has issues with modelling contact discontinuities, due to the surface tension issue I mentioned. These are potentially important in cosmological simulations, where you get large density contrasts, for example. The classical test seems to be the (rather photogenic) Kelvin-Helmholtz instability, where two initially laminar flows in different fluids are set up to slip past one another. What should happen is that you get turbulent waves forming at the interface of the fluids, but the numerical “tension” in SPH models stops this from happening. Daniel Price took up the gauntlet and attempted to show that SPH can be fixed by properly taking into account discontinuities in the density. First of all, he showed that a standard SPH code will happily produce realistic Kelvin-Helmholtz instability if the two fluids have comparable densities – the problem is not that SPH can’t deal with instability, it’s that it can’t deal with density discontinuities, since SPH normally assumes that the density is differentiable at boundaries. By switching a viscosity term on and off depending on the situation (e.g. Price 2011; Read and Heyfield 2011), Daniel claims that SPH can be made viable for cosmological simulations again. Volker disagreed in the questions afterwards, and I think they ended up “agreeing to disagree”.
There were also talks by Romain Teyssier and Brian O’Shea discussing the differences between different types of simulation, and comparing different codes (amongst other things). There have been a number of code comparisons in recent years, for a number of different physical situations (e.g. Frenk et al. 1999; Agertz et al. 2007; Scannapieco et al. 2011), which seem to show that most codes get broadly similar results on larger scales (for massive dark matter halos), but have divergent results at smaller scales, where baryonic physics becomes important.
Broadly similar sentiments were expressed by Rob Thacker, who talked about modelling AGN feedback. To do feedback, you need an extremely large dynamical range – in feedback studies, processes from galaxy scales down to the supermassive black hole event horizon (or thereabouts) are important. Adaptive Mesh Refinement codes (like Teyssier’s RAMSES) can get something like the dynamical range required, and not-bad results on large scales, but it sounds like there are lots of problems and uncertainties with the models of the small-scale physics. (Fun RAMSES fact: It uses Peano-Hilbert space-filling curves to do “domain decomposition” – deciding which particles are governed by which processor for massively parallel execution.) He also noted, somewhat surprisingly, that black hole advection might be important during galaxy mergers (the BH gets moved around under the gravitational pull of nearby matter, thus experiencing a changing local density, which has an effect on the accretion rate, and so on), even though the BH is extremely massive compared to anything nearby. A final interesting nugget came from the audience – apparently, during BH mergers you should expect to lose about 5% or so of the total BH-BH system mass in the form of gravitational waves, but Rob was confident in saying that current AGN simulations are nowhere near having to worry about 5% effects like that!
Lots of other interesting things went on today, but I’ll wrap up with a couple of brief notes. Ilian Iliev talked about simulating reionisation, and the challenges involved in getting photoionisation right. Apparently, many people use a purely local form for the photoionisation rate equation, which does not take into account finite volume effects, and thus does not properly conserve photon number. He wrote down the correct version, which includes an extra term in the frequency integral. Romain also mentioned that many ray-tracing codes assume that the speed of light is infinite, which is often a good approximation, but which sometimes isn’t, and gets forgotten. This can lead to up to a factor of two discrepancy in some calculations. Finally, Brian mentioned a paper by Tasker et al. (2008), which lists a suite of useful test cases for numerical simulation codes. Many of these tests were wheeled out today.
Andreas Bauer, & Volker Springel (2012). Shocking results without shocks: Subsonic turbulence in smoothed particle hydrodynamics and moving-mesh simulations MNRAS DOI: http://arxiv.org/abs/1109.4413v1