Overview -- Numerical Relativity Using CUDA C/C++ Is Easier Than You Think!

Simulation results of a binary black hole head-on collision run on a GPU based home gaming computer

Not too long ago, the only way to numerically solve the Einstein-Hilbert field equations required the use of a supercomputer that ran a hundred or more CPUs in parallel for days at a time. In 2005, Frans Pretorius used a cluster of 125 CPUs at the University of Texas to run the first every numerical simulation of the merger of an orbiting binary pair of black holes. His simulation took 30 days to complete.

Univ. of Texas Advanced Computer Center - the Stampede Cluster

link to Frans Pretorius Home Page

For massive simulations, like cosmology models, large clusters and packaged routines are still the way to go. Bu, for smaller simulations, like head-on binary black hole simulations, this is no longer the case.

Since the wide-spread use the Graphics Processor Units (GPU) in home computers, the average Joe has access to systems that run at TFLOP (Trillions of Floating-Point Operations per second) rates.

HP Omen Home Gaming Computer w/ 11-TFLOPS rated NVIDIA GeForce 1080-Ti

NVIDIA GeForce 1080Ti GPU

The problem is that many physicists I have spoken with who do numerical relativity for a living are not at all familiar with programming CUDA C/C++. This is a real bummer, because it means one is limiting oneself to either using canned routines or only clusters that use CPUs.

Every serious researcher interested in the best simulation is going to want to go beyond canned routines. Better yet, I  would hope that those  folks who maintain canned numerical relativity routines would be interested in porting their code to CUDA to exploit the massive speed ups one can get on GPU clusters.

And so, I've developed this tutorial to help people who want to utilize the power of GPUs in their own work.

My view is, if I can do it, certainly the people who have made a career from numerical relativity can learn to program in CUDA C/C++ as well.

Methods abound, but for the purposes of this example, I will focus on use of the BSSN moving puncture method with dynamic gauge variables in Cartesian coordinates using a 3D second order finite difference method for spatial derivatives and a Euler finite difference for evolution in time. 

The simulation results presented in this video were completed on a HP Omen gaming computer that is equipped with a NVIDIA GeForce 1080-Ti GPU card, which is rated at 11 TFLOPS. The language used in the simulation is CUDA C/C++. The total simulation time was about 8 hours or so. And, the code is hardly optimized.

Yes, I chose simplicity for the purpose of this example. The point is that if one were to learn how to program a simple example, then hopefully one can take those learnings and apply them to programming large clusters of GPUs.

No doubt these methods can be improved, and it is the hope of the author that the reader will do so in his or her own version. However, the purpose here is to demonstrate how to utilize the power of the GPU to perform a realistic numerical relativity example. In that sense, the simpler the method the better. In this case, we will examine how to use CUDA C/C++ to simulate a head-on collision of two black holes.

The video below is the result of one such simulation.Again, it took about 8-hours to run to completion. The graphs are of the W variable used in the moving puncture BSSN formulation, which can be thought of as the strength of the gravitational field.

References

[1] Shibata, Masaru. Numerical Relativity :100 Years of General Relativity. World Scientific Publishing Company (November 5, 2015).

 [2] CUDA C Programming Guide, https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html

[3] Duane Stori and Mete Yurtoglu. CUDA for Engineers: An Introduction to High-Performance Parallel Computing. Addison-Wesley Professional; 1st edition (November 12, 2015).

Next: Background Material on General Relativity -->

Table of Contents
Overview
Background Material on General Relativity
Cracking the Code
Moving Puncture BSSN Evolution Equations
Dynamical Gauge Equations
Initial Conditions for Head-On Collision of Two Black Holes
Kreiss-Oliger Spatial Filtering
Structures Used to Implement on GPU Machine
Implementation on CPU example
Implementation on GPU Example
Finite Differences on a GPU Machine
Variable Glossary

Background Material on General Relativity

Numerical Relativity is the use of a computer to simulate the Einstein-Hilbert field equations. These relate the curvature of space-time, denoted by the trace-free Ricci curvature tensor, and the cosmological constant to the stress-energy-momentum tensor:

It also appears with the Einstein tensor (which is the same thing as the trace-free Ricci curvature tensor)

These equations can be categorized as nonlinear wave equations in 4-dimensional space-time. 

Unlike the space-time of Newton, or even of special relativity, the space-time of general relativity can be curved. This curved space-time has the general form of a pseudo-Riemman manifold, that has the general metric:

Almost immediately after the discovery of the field equations of gravity, Schwarzchild published his solution to a radially symmetric space-time, the Schwarzchild metric (in units in which the speed of light, c, is equal to unity):

Several other analytical solutions soon followed, including Einstien’s own steady-state cosmological model, de Sitter’s exponentially growing vacuum solution, and Friedmann’s more general cosmological model (also called the Friedmann–Lemaître–Robertson–Walker cosmological model.) Except for Einstein’s steady-state universe model, versions of these models are still used as approximations to the cosmology of our universe today.

Unfortunately, to go beyond symmetric geometries to more realistic systems, the Einstein-Hilbert field equations need to be solved numerically.

Next: Cracking the Code -->

< -- Back: Overview

Cracking the Code

Well, if the goal is to simulate the merger of two black holes, then indeed there is one very big simplification that can be made. Black holes, though they have mass, are actually considered to be vacuum solutions. That is, the stress-energy-momentum tensor (the right-hand side of the Einstein-Hilbert field equations) is zero.

Another simplification can be made by noting that the cosmological constant, Λ_μν is quite small. In SI units, it is measured to be on the order of 10^-52m^-2.Though it might contain the majority of the energy in the observable universe, this is only when taken on vast, cosmological scales of the Mega-parsec range. In the few hundred billion-km or so that would be used to simulate even the largest black hole mergers, the effect of dark energy is thought to be minuscule by comparison to the other terms. Thus, for our purposes, setting Λ_μν to zero is a safe bet.

Thus, the version of the Einstein-Hilbert field equations that we would to wrestle with is:

Okay, this is beginning to look easy!

Relatively speaking, yes. Still, let's take a look at what the vacuum version is made up of before jumping to any such conclusion. First, note that the tensor equations we are left with have two indices. And, for four-dimensional space-time, that means there are 4 x 4 = 16 interdependent equations to solve.

Well, we can do a little better than that, since the Einstein tensor, E_μν , is symmetric (E_μν = E_νμ ) there are then just 10 non-redundant nonlinear partial differential equations to deal with. This is also true of the Ricci curvature tensor, Rμν .

Then, the system of equations we would like to solve looks like:

This looks like progress. But, just what is the Ricci curvature tensor? Ricci curvature is a measure of how much the local curvature deviates from Euclidean (i.e., non-curved) space in a Riemannian manifold. It can be calculated as a twice contraction of the Riemannian curvature, Rμbνc, with the metric g^bc

In this last equation, the Einstein summation convention is used, which is to say it can also be written in perhaps a more standard mathematical form as:

So then, what is the Riemannian curvature? The Riemannian curvature is a measure of curvature for manifolds with dimensions greater than 2 and is calculated as (note that α is used to indicate summation in the last two terms on the right hand side of the equation below):

where Γ ^λ _μσ  is the Christoffel symbol, which is calculated as:

So, when you unravel the Einstein-Hilbert field equations into their raw form, you find that they represent the relation between the space-time metric (g_αβ or g^αβ) and the slope of the space-time metric (∂_βg_μν.) So, given the value of the metric at all points in our 4-dimensional grid at some initial time, we should be able to figure how, based on the slope of the metric, it evolves.

(Actually, if you had the metric (and its slope) for all points in a 4-dimensional grid, you'd have the answer for all time. Which, I have to admit, sounds handy.)

It turns out that is a little tricky. For example, what does "initial time" in a 4-dimensional space-time even mean? 

Maybe it would be nice to formulate the general relativity field equations into a more standard initial value problem. 

This is the goal of a method called the 3+1 numerical relativity technique, which separates out time from space time (hence "3" spatial dimensions "+1" time dimension.)

Next: Moving Puncture BSSN Evolution Equations -->

<-- Back: Background on General Relativity
<--Back: Overview

ADM 3+1 Formalism

The story goes, Arnowitt, Deser and Misner weren't necessarily looking for a numerical relativity method when they developed their Hamiltonian version of the Einstein-Hilbert field equations. No, they were looking for a way to develop a quantum theory of gravity.

But, their clever method of treating 3-dimensional foliations of space within the 4-dimensional space-time manifold allowed those interested in numerical relativity to format the general relativity field equations as an initial value problem.

Foliation here is meant to be the breaking up 4-dimensional space-time into non-intersection 3-dimensional spatial hypersurfaces. (Think of the 2-dimensional leaves of a tree forming separate surfaces in 3-dimensional space.) The drawing below shows two nearby hypersurfaces a "distance" Δt away from each other in time. The lower one is denoted as Σ_Δt and the upper one as Σ_t+Δt

Note that the vector that represents the time dimension t above in the drawing (the vector from point 'O' to point 'B') may actually not be parallel to the normal (the vector from point 'O' to point 'A') to the surface Σ_Δt. This means that in some increment of time, t , the distance between the two hypersurfaces along the normal to Σ_Δt  is actually α times Δ t instead of t. The vector between where the normal lands on Σ_t+Δt  and Σ_t+Δt  is  β.

Another name for α is the lapse function (as in a lapse of time) and another name for the offset vector β between point 'A' and point 'B' is the shift function. Collectively, α and β are also known as the gauge variables.

BSSN Moving Puncture Evolution Equations

There are several great texts on numerical relativity. The one I used primarily is M. Shibata’s book [1]. His work lays out in great detail many different aspects of numerical relativity in general, including the BSSN moving puncture method in Cartesian coordinates. 

The basic evolution equations (in vacuum) are:

where:

[1] Shibata, Masaru. Numerical Relativity (100 Years of General Relativity). World Scientific Publishing Co. Pte Ltd.

Next to Dynamical Gauge Variables -->

<--Back: Cracking the Code

<-- Back to Overview

Dynamical Gauge Equations

This version uses dynamical gauge variables, so there are several more equations to code as well.

Initial Conditions for Head-On Collision of Two Non-spinning Black Holes

Additionally, for this example, since this is a head-on collision of two non-spinning black holes, our initial conformal parameter will be calculated as:

Where  m_i is the mass associated with each black hole, is the initial location of each black hole and the location of  each point is

Kreiss-Oliger Spatial Filtering

Because of high frequency numerical noise, this version also includes Kriess-Oliger spatial filtering. This filter is implemented as:

Next: Initial Conditions -->
<-- Back to BSSN Moving Puncture Evolution Equations
<-- Back to Overview

Initial Conditions for Head-on Collision of Two Non-Spinning Black Holes

As in the rest of this example, we choose a simple case, that of two black holes, neither of which has spin or boost (angular or linear momentum.)  For the moving puncture version of BSSN, the Brill-Lindquist initial data method is used.

In this example, since this is a head-on collision of two non-spinning black holes, our initial conformal parameter:

will be calculated (again, using the Brill-Lindquist method) as:

                                                                                                                                

Where:

is the mass associated with each black hole, 

is the initial location of each black hole 

is the location of  each point


Next: Kreiss-Oliger Spatial Filtering -->

<-- Back: Dynamical Gauge Variables

<--Back: Overview

Kreiss-Oliger Spatial Filtering

Because of high frequency numerical noise, this version also includes Kriess-Oliger spatial filtering. This filter is implemented as:

Note that, since we are using  a 3-dimensional spatial grid, this filter is taken in all three spatial dimensions. That is to say, we perform the filter, for example, first on the x-dimension, then in the y-dimension, and then in the z-dimension. (Obviously, the order is not important here.)

Next: Structures Used to Implement on GPU -->

<-- Back: Initial Conditions

<-- Back: Overview