What is a numerical simulation?
In astrophysics most processes can be described by a set of equations. A simple example are the equations of motion, which describe how a body moves under the effect of forces that produce acceleration. Not all equations are easy to solve analytically, and some of them don’t even have analytical solutions, in which case we have to solve them numerically. To add more complications sometimes there are so many processes involved that is really hard to handle all equations. Here is where we appeal to numerical simulations. In simple words, a numerical simulation is a computational program that given an initial condition is able to solve the equations governing the problem and calculate the evolution of the system.
Why use simulations?
In principle the definition of numerical simulations may sound not too difficult, but in reality a system is determined by many equations that need to be solved at the same time. Things get more complicated if we want a better resolution and we use too many particles or too many grids. In practice super computers and multiple processors are needed to do a detailed simulation. The power of simulations resides on that they allow us to study processes that we are unable to observe. Reasons for that are many. For example some processes such as collision of galaxies take a long time and we are not able to watch the whole process with telescopes. Another interesting example is the formation of the first structures in the Universe. A very famous simulation of large scale structure is the Millenium simulation. In the context of shocks, give us insights on how shocks are formed, for example through AGN feedback or Supernovae explosions
There are two main techniques to develop a code:
- Adaptive Mesh Refinement (AMR): Based in an Eulerian hydrodynamics method for solving equations, it consists of fixed grids sampling the spatial domain of the simulation. The properties of the fluid are stored in each grid according to their position in space. Grids can be subdivided in subgrids in order to get better resolution. This refinement technique can depend on various factors (e.g. reaching a density threshold) but the most important is the Truelove condition which says that the Jeans mass has to be resolved by at least four cells to avoid numerical fragmentation. In AMR codes the Euler equations are usually solved using finite-difference methods, but recent codes also include finite-volume and finite-element methods. Examples of AMR codes are Enzo (webpage,paper), Flash (webpage,paper), Ramses (webpage,paper), Orion (webpage), Athena (webpage,paper), Castro (webpage,paper)
- Smoothed Particle Hydrodynamics (SPH): Based in an Lagrangian hydrodynamics method for solving equations, it consists in particles sampling the fluid. That means that particles move with the fluid, and hence their Lagrangian nature. Since there is only a limited number of particles the whole main is not covered. To extrapolate the fluid properties to regions that are not sampled by the particles a smoothing kernel is used. The resolution depends on the mass of each particles, the smaller the mass the better the resolution. Similar to the Truelove condition, in SPH the mass resolution has to be less than the Jeans mass to avoid numerical fragmentation. In SPH the fluid equations are solved based on the Euler-Lagrange equations . Examples of SPH codes are Gadget (webpage,paper), Gasoline (paper), Hydra (paper)
Recently a third method has been developed called Moving Mesh, which is a mixture of AMR and SPH. We will not give details about the code since it is still private, but if you are curious you can read the paper that describes the code.