by Devon Powell

A complicated calculation

In the early summer of 1945, physicist Bob Christy asked fellow physicist Richard Feynman to carry out a task as quickly as possible. The deadline was the Trinity nuclear test, the first nuclear bomb and the culmination of years of secret work by Manhattan Project scientists. The task was to predict the total energy that would be released by the Gadget device, the prototype implosion bomb designed at Los Alamos.

Although this type of calculation had been done before, the rush caused Feynman to drop everything and rally his troops. In those days, mathematical operations were carried out by literal human “computers” either working by hand (i.e. looking at a table of logarithms) or using mechanical arithmetic machines. This meant that each human computer could perform one operation every few seconds. Time management was of the utmost importance.

The resulting system was, in typical Feynman style, elegant. Predicting the bomb yield involved many high-level calculations, each consisting of many computer operations. These calculations did not all depend on each other, so some could be worked on independently. Results of individual operations were passed from human computer to human computer on cards, which were color-coded by the calculation they belonged to. In addition to the colors allowing operations to be done in parallel, they also allowed cards from different calculations to be mixed. This meant that a queue of cards from any number of different calculations could be fed through the same specialized mathematical operation, and sorted back out afterwards. Jobs of different priorities could also be reordered. Errors could also be quickly identified and traced back through a stack of cards, so that downtime for error correction was minimized. This system today might be referred to as an algorithm.

(Note: important technical terms are in italics, and can be found in a short glossary at the end of this post.)

Feynman and his team finished the job. The types of physics they needed to simulate the bomb are typical of many science applications in astrophysics and cosmology. Some of the computational techniques they pioneered to help simulate them are still in use at the world’s largest supercomputers. Hydrodynamics (the study of moving fluids), radiation transport (modeling the propagation of light or radiation through a fluid), and nuclear physics (modeling interactions between radiation and atomic nuclei) are all fields of active research in modern astrophysics. Together, hydrodynamics, radiation transport, nuclear physics, chemistry, and gravity are important for studying nuclear fusion in stars and supernovae, interstellar gas, ionization, accretion around black holes, and much more.

What does it mean to simulate a physics problem on a computer?

Astrophysical processes are described by partial differential equations (PDEs). This is a special type of equation describing how physical systems change over time. Here are several examples:

Don’t worry about understanding these equations if you are unfamiliar with the pieces of them—I’m just showing you to give you a sense of their complexity (the actual solutions would occupy many more pages, yet). They are called non-linear PDEs by physicists, and they are very difficult to solve by hand. Analytic solutions to these equations (mathematical expressions that can be written down succinctly) are rare and are often uninteresting for learning about the complex physics in space. This means that in most cases, the only way to study a complicated nonlinear system in detail is to simulate it on a computer.

In fact, solving nonlinear PDEs is so difficult that when we simulate them on the computer, we are actually solving for an approximation to the true solution. This is by necessity, since computers only have a finite amount of memory and processing power. So, we reduce the amount of information the computer needs to hold by discretizing the data. This means that we break the simulation into a large, but manageable, set of numbers in computer memory. The way in which this is done must be derived mathematically so that the discrete version of the simulation is consistent with the original equations. There are many ways to do this for different equations, but the important thing to remember is that we can run a simulation on some finite data set meant to approximate the true solution.

Clearly, though, an approximation is not a true solution to a problem. This is a philosophical challenge often posed to computational physicists: Why bother simulating something that is just an approximation? As it turns out, the Lax Equivalence Theorem states that if we did have an infinitely large computer, and we discretized our simulation into an infinitely large memory bank, then we would in fact get the exact solution. This is true as long as the discretized simulation is mathematically consistent with the original physics PDE. Furthermore, the more information we store on a real computer with finite memory, the closer to the true solution our approximation becomes. This is a very powerful phenomenon called convergence. In a nutshell, it means that a higher fidelity simulation will produce a result closer to the true solution.

Imagine two digital photos side-by-side: the one taken with a 10-megapixel camera will look much better than the one taken with a first-generation iPhone. It’s not as perfect as looking at the actual scene with your eyes, but it’s extremely close. In the same way, as we apply more and more memory and processing power to our simulation (of, say, the Milky Way, as an example) we will continue to develop a much clearer scientific picture.

What I use this for

My research interests lie in using computational geometry to create more physically consistent discretizations for astrophysical systems by taking advantage of 3D solid geometry calculations. During my graduate career, I have written or improved upon computer codes that solve all three of the above PDEs: the Euler equations for hydrodynamics, the Vlasov-Poisson equation for gravity, and the radiative transfer equation. These are three very different types of physics, so different algorithms are needed to solve them. The latter two are most relevant to my work, so I will focus on them.

The Vlasov-Poisson equation describes how dust or dark matter (any pressureless fluid) moves under the influence of its own gravity. This is the physical system that N-body simulations approximate in order to study the distribution of matter on cosmic scales. I published a paper with my advisor, Tom Abel, on the use of a new dark matter discretization to create much smoother and clearer pictures of the matter density distribution in the Universe. This was the topic of a previous KIPAC blogpost, "Where the invisible things are: Probing the phase space structure of dark matter". Our method relies technically on being able to compute the exact volume contained in the intersection between a tetrahedron and a cube. Researching this seemingly simple problem exposed me to a world of computational geometry that had never been applied to physics calculations.

Below is a figure showing the time evolution of the cosmic density field using this new discretization. The simulation is smoother and more physically consistent with the original Vlasov-Poisson equation than traditional N-body methods. However, the new method requires many more individual computer operations.

I have most recently been applying the same computational geometry algorithm used for the Vlasov-Poisson system to the radiation transport problem. The radiative transfer equation is very interesting in astrophysics and nuclear physics. The first stars in the Universe drastically altered the chemistry of the cosmos with their ultraviolet radiation, a process called “reionization” that we would like to understand better. The radiative transfer equation above looks complicated, but has a rather straightforward physical meaning: radiation in a vacuum should propagate freely, while radiation passing through matter should be absorbed or scattered. The details of how this new method works are largely based on ray-tracing methods. Our improvement on ray tracing was to pass radiation through precisely calculated intersection volumes between grid cells and a beam of radiation. This gives a much more accurate calculation of the amount of radiation that should be absorbed or allowed to pass through the matter in a grid cell. There are many types of radiation that can be simulated using this method. Here, “radiation” refers the the ultraviolet photons (UV light) emitted by the first stars and galaxies. However, radiation transport is also used to track nuclear radiation (neutrons, protons, and gamma-ray photons) inside for example, supernovae and nuclear bombs.

We again found that by using the more physically consistent discretization given by the geometry of the problem, we were able to get very accurate results while using less memory than previous radiation transport methods. Below is a picture from a simulation that used ~200 data structures to describe the radiation field, rather than many thousands for the previous ray-tracing methods.

Both of these projects improve upon the accuracy of simulations by making the numerical algorithms more consistent with the equations describing them. For the Vlasov-Poisson system, we create accurate density maps by finding the exact intersections between tetrahedral mass elements and cubical grid cells. For radiation transport, we pass light through matter by calculating the geometrically exact overlap between the grid cells and the radiation beam. These methods improve on N-body codes and ray-tracing codes, respectively, by more accurately sampling the simulation space using 3D solid geometry rather than point particles or 1D rays. This method is therefore more consistent with the original physical system we want to study. This consistency comes at the cost of computer operations, since the geometric calculations involved are relatively complex. However, supercomputers have increased in power by many orders of magnitude since the first ones were built. This is a result of more and faster processors as well as new options for accelerating code on graphics processing units (GPUs). More and more available computing power will enable more accurate but more complex calculations like the ones I described here to be run in reasonable amounts of time for large-scale simulations.

We have come a loooong way, now.

In our Manhattan Project anecdote, we noted that each human could perform one operation every few seconds. We measure the speed of a computer in FLoating point Operations Per Second, or “flops”. So, a human computer runs at less than a flop. Your laptop can do a few billion (few x 10^9) flops. The “fastest” supercomputers are pushing the petaflop to exaflop range (10^15 to 10^18), though supercomputers are really just a bunch of regular computers linked together in a cooperative network.

Compared to modern astrophysics calculations, Richard Feynman’s bomb yield was a relatively simple one. Simulation techniques have continued to develop rapidly since then. New algorithms have greatly increased our ability to run fully three-dimensional, multi-physics simulations with higher and higher resolution. This has given us a much better understanding of astrophysical processes, and will continue to do so as supercomputers become faster and faster.

Despite the huge increase in speed, supercomputers are fundamentally built to do the same thing as the Los Alamos team. They have memory in which to order and store results, they have processors on which to execute different mathematical operations, and they have network links that allow them to share results between processors. The difference is that modern scientists tell the computer what to do by writing software, rather than tediously coordinating the movement of colored cards by hand.

And we are very, very thankful for that!

Glossary

Algorithm: A set of instructions given to the computer in order to execute a calculation. Algorithms can be simple or complex depending on the application.

Approximation: A solution that is close, but not exactly equal to, the exact solution of an equation. Physics computations on a computer all solve for an approximation to the true solution, but better and faster computers can get ever closer to the true solution. Consistent: When deriving the discrete form of a physics equation, one must make sure that the derivation is correct with regard to the original equation, and will give an approximation close to the exact solution.

Convergence: If we solve for an approximation to the discrete version of a non-linear PDE, then the approximation will be closer and closer to the true solution as we use more memory and processing power. This is true if the discrete equation is consistent with the original.

Discretize/discrete: To break up the information describing a computer simulation (i.e. the density field) into a large but finite set of numbers. The modified physics equations used to solve this finite, approximate system are known as the discrete form of the equations.

Non-linear PDE: Stands for “non-linear partial differential equation.” This is a special type of physics equations that describes how complicated systems change over time.

References and further reading

• The Pleasure of Finding Things Out, by Richard Feynman.
• Numerical methods for ordinary differential equations
• Lax equivalence theorem
• Stability in numerical differential equations
• Convergence, Consistency, and Stability (technical PDF)
• Properties of numerical methods (technical PDF)