Ray Tracing in Curved Space Time
An example of a HPC task is to calculate the trajectories of photons near a black hole. The extreme physics means that they no longer travel in straight lines due to the curvature of space-time. So object near a black hole will be gravitationally lensed, when viewed from afar. Here, we will investigate how to calculate the observed shape of a thin disk surrounding the hole.
Having a thin disk around a black hole isn't too outlandish. When gas orbits around a black hole, the conservation of angular momentum means that friction forces will eventually convert it into an "accretion disk". In common situations, that disk will be thin compared to its extent. Astronomers using telescopes to observe such systems need to understand their results, so having a theoretical model of what the simplest possible case will produce is quite helpful.
The first step is to model the rays of light as they travel. The problem here is that gravity curves space-time near a black hole. Using simple geometry to calculate intersections no longer is possible. (You might be able to describe the trajectories in terms of elliptic functions... but that doesn't help all that much for more complex problems.) Instead, we can describe the motion of the photon as a set of coupled differential equations in Boyer-Lindquist Coordinates:
/* Coupled differential equations describing motion of photon */
static void geodesic(double *y, double *dydx)
{
double r, theta, pr, ptheta;
r = y[0];
theta = y[1];
pr = y[4];
ptheta = y[5];
double r2 = r*r;
double twor = 2.0*r;
double sintheta, costheta;
sincos(theta, &sintheta, &costheta);
double cos2 = costheta*costheta;
double sin2 = sintheta*sintheta;
double sigma = r2+a2*cos2;
double delta = r2-twor+a2;
double sd = sigma*delta;
double siginv = 1.0/sigma;
double bot = 1.0/sd;
/* Prevent problems with the axis */
if (sintheta < 1e-8)
{
sintheta = 1e-8;
sin2 = 1e-16;
}
dydx[0] = -pr*delta*siginv;
dydx[1] = -ptheta*siginv;
dydx[2] = -(twor*a+(sigma-twor)*L/sin2)*bot;
dydx[3] = -(1.0+(twor*(r2+a2)-twor*a*L)*bot);
dydx[4] = -(((r-1.0)*(-kappa)+twor*(r2+a2)-2.0*a*L)*bot-2.0*pr*pr*(r-1.0)*siginv);
dydx[5] = -sintheta*costheta*(L*L/(sin2*sin2)-a2)*siginv;
}
Where we have used some constants of the motion to simplify things a bit: L, for the angular momentum in the phi direction. kappa, for the something related to "Carters Constant", which in turn is related to the angular momentum in the theta direction. And finally we have set the energy at infinity of the photon to unity to normalize. The resulting math isn't too complex. The only thing to worry about is if a photon gets too close to the rotational axis of the black hole. If it does, our coordinate system becomes degenerate, and we need to worry a bit about divisions by zero. (The same problem happens when you are at the north pole and want to know which longitude you are at...)
To start the photons off, we need some initial conditions:
/* Initial Conditions for Ray */
static void initial(double *y0, double *ydot0, double x, double y)
{
y0[0] = r0;
y0[1] = theta0;
y0[2] = 0;
y0[3] = 0;
y0[4] = cos(y)*cos(x);
y0[5] = sin(y)/r0;
double sintheta, costheta;
sincos(theta0, &sintheta, &costheta);
double cos2 = costheta*costheta;
double sin2 = sintheta*sintheta;
double rdot0 = y0[4];
double thetadot0 = y0[5];
double r2 = r0 * r0;
double sigma = r2 + a2*cos2;
double delta = r2 - 2.0 * r0 + a2;
double s1 = sigma - 2.0 * r0;
y0[4]= rdot0*sigma/delta;
y0[5]= thetadot0*sigma;
ydot0[0] = rdot0;
ydot0[1] = thetadot0;
ydot0[2] = cos(y)*sin(x)/(r0*sin(theta0));
double phidot0 = ydot0[2];
double energy2 = s1*(rdot0*rdot0/delta+thetadot0*thetadot0)
+ delta*sin2*phidot0*phidot0;
double energy = sqrt(energy2);
/* Rescale */
y0[4] = y0[4]/energy;
y0[5] = y0[5]/energy;
/* Angular Momentum with E = 1 */
L = ((sigma*delta*phidot0-2.0*a*r0*energy)*sin2/s1)/energy;
kappa = y0[5]*y0[5]+a2*sin2+L*L/sin2;
/* Hack - make sure everything is normalized correctly by a call to geodesic */
geodesic(y0, ydot0);
}
Where we have an observer at some distant location described by r0,theta0 looking back at the disk. The photons are parameterized by x, y in the "view port" of the observer. By scanning over x and y, we can build up an image once we know where each photon goes.
Now that we have the math out of the way, the next step is to worry about calculating the motion. Since we have a set of coupled ordinary differential equations, we use a Runge-Kutta algorithm for this. It numerically solves the equations, giving us the position and direction of the photon at discrete points along its path. The R-K algorithm will dynamically tune the sizes of the steps it takes in order to maintain error within some given tolerance. We take a standard R-K routine, and rewrite it so that the compiler can easily vectorize the loops:
static void rkstep(double *y, double *dydx, double h, double *yout, double *yerr)
{
int i;
double ak[N];
double ytemp1[N], ytemp2[N], ytemp3[N], ytemp4[N], ytemp5[N];
for (i = 0; i < N; i++)
{
double hdydx = h * dydx[i];
double yi = y[i];
ytemp1[i] = yi + 0.2 * hdydx;
ytemp2[i] = yi + (3.0/40.0) * hdydx;
ytemp3[i] = yi + 0.3 * hdydx;
ytemp4[i] = yi -(11.0/54.0) * hdydx;
ytemp5[i] = yi + (1631.0/55296.0) * hdydx;
yout[i] = yi + (37.0/378.0) * hdydx;
yerr[i] = ((37.0/378.0)-(2825.0/27648.0)) * hdydx;
}
geodesic(ytemp1, ak);
for (i = 0; i < N; i++)
{
double yt = h * ak[i];
ytemp2[i] += (9.0/40.0) * yt;
ytemp3[i] -= 0.9 * yt;
ytemp4[i] += 2.5 * yt;
ytemp5[i] += (175.0/512.0) * yt;
}
geodesic(ytemp2, ak);
for (i = 0; i < N; i++)
{
double yt = h * ak[i];
ytemp3[i] += 1.2 * yt;
ytemp4[i] -= (70.0/27.0) * yt;
ytemp5[i] += (575.0/13824.0) * yt;
yout[i] += (250.0/621.0) * yt;
yerr[i] += ((250.0/621.0)-(18575.0/48384.0)) * yt;
}
geodesic(ytemp3, ak);
for (i = 0; i < N; i++)
{
double yt = h * ak[i];
ytemp4[i] += (35.0/27.0) * yt;
ytemp5[i] += (44275.0/110592.0) * yt;
yout[i] += (125.0/594.0) * yt;
yerr[i] += ((125.0/594.0)-(13525.0/55296.0)) * yt;
}
geodesic(ytemp4, ak);
for (i = 0; i < N; i++)
{
double yt = h * ak[i];
ytemp5[i] += (253.0/4096.0) * yt;
yerr[i] -= (277.0/14336.0) * yt;
}
geodesic(ytemp5, ak);
for (i = 0; i < N; i++)
{
double yt = h * ak[i];
yout[i] += (512.0/1771.0) * yt;
yerr[i] += ((512.0/1771.0)-0.25) * yt;
}
}
The driver routine that determines the step sizes can be more standard:
static double rkqs(double *y, double *dydx, double htry, double escal, double *yscal, double *hdid)
{
int i;
double hnext;
double errmax, h = htry, htemp;
double yerr[N], ytemp[N];
while (1)
{
rkstep(y, dydx, h, ytemp, yerr);
errmax = 0.0;
for (i = 0; i < N; i++)
{
double temp = fabs(yerr[i]/yscal[i]);
if (temp > errmax) errmax = temp;
}
errmax *= escal;
if (errmax <= 1.0) break;
htemp = 0.9 * h / sqrt(sqrt(errmax));
h *= 0.1;
if (h >= 0.0)
{
if (htemp > h) h = htemp;
}
else
{
if (htemp < h) h = htemp;
}
}
if (errmax > 1.89e-4)
{
hnext = 0.9 * h * pow(errmax, -0.2);
}
else
{
hnext = 5.0 * h;
}
*hdid = h;
memcpy(y, ytemp, N * sizeof(double));
return hnext;
}
This allows us to follow the photon where-ever it goes. The next trick is to worry about imaging the disk. The problem here is that the R-K routine will give us discrete steps as output. These steps most likely will not ever exactly hit the disk, and instead will skip over it. What we need is some way to tune the size of the step in order to hit the disk to within some tolerance. Since space-time is curved, there is no easy way of guessing the size of step needed. Instead, we use a binary search technique by noticing which hemisphere the photon happens to lie in.
static void binarysearch(double *y, double *dydx, double hbig)
{
double hsmall = 0.0;
int side;
if (y[1] > M_PI/2.0)
{
side = 1;
}
else if (y[1] < M_PI/2.0)
{
side = -1;
}
else
{
/* Already at the equator */
return;
}
geodesic(y,dydx);
while ((y[0] > Rhor) && (y[0] < r0) && (side != 0))
{
double yout[N], yerr[N];
double hdiff = hbig - hsmall;
if (hdiff < 1e-7)
{
rkstep(y, dydx, hbig, yout, yerr);
memcpy(y, yout, N * sizeof(double));
return;
}
double hmid = (hbig + hsmall) / 2;
rkstep(y, dydx, hmid, yout, yerr);
if (side * (yout[1] - M_PI/2.0) > 0)
{
hsmall = hmid;
}
else
{
hbig = hmid;
}
}
}
With that done, it is an easy matter to fire rays from our observer, and then see if they hit the disk or not. In order to image the disk, we break it up into a checker-board pattern. The following routine then will return the rgb value of the pixel corresponding to a ray. If the ray avoids the disk, and either escapes to infinity, or hits the black hole, we colour the pixel black. Otherwise, we will colour it a blue shade based on where in the disk it hit.
static void fire_ray(unsigned char *rgb, int x1, int y1)
{
double htry = 0.5, escal = 1e11, hdid = 0.0, hnext = 0.0;
double range = 0.0025 * Rdisk / (size - 1.0);
double y[N], dydx[N], yscal[N], ylaststep[N];
int side;
int i;
initial(y, dydx, (x1 - (size + 1.0) / 2) * range, (y1 - (size + 1.0) / 2) * range);
while (1)
{
memcpy(ylaststep, y, N * sizeof(double));
geodesic(y, dydx);
for (i = 0; i < N; i++)
{
yscal[i] = fabs(y[i]) + fabs(dydx[i] * htry) + 1.0e-3;
}
if (y[1] > M_PI/2) side = 1;
else if (y[1] < M_PI/2) side = -1;
else side = 0;
hnext = rkqs(y, dydx, htry, escal, yscal, &hdid);
if ((y[1]-M_PI/2)*side < 0)
{
memcpy(y, ylaststep, N * sizeof(double));
binarysearch(y, dydx, hdid);
/* Did we hit the disk? */
if ((y[0] <= Rdisk) && (y[0] >= Rmstable))
{
unsigned p1 = 4.0 * (y[0] - Rmstable) / (Rdisk - Rmstable);
unsigned p2 = floor(y[2] * 6.0 / M_PI);
if ((p1 ^ p2) & 1)
{
rgb[0] = 255;
rgb[1] = 128;
rgb[2] = 128;
}
else
{
rgb[0] = 255;
rgb[1] = 0;
rgb[2] = 0;
}
return;
}
}
/* Inside the hole, or escaped to infinity */
if ((y[0] < Rhor) || (y[0] > r0))
{
rgb[0] = 0;
rgb[1] = 0;
rgb[2] = 0;
return;
}
htry = hnext;
}
}
So now, we just need to worry about outputting the rgb image data. The easiest way of doing that is to use the Targa (tga) format. It has a very simple header, and raw rgb follows it. A couple of routines to output the header are:
/* Write all the data, or fail and exit */
static void full_write(int fd, void *buf, size_t size)
{
while (size)
{
ssize_t did = write(fd, buf, size);
if (did == -1)
{
if (errno == EINTR) continue;
errx(1, "Error writing file: %s\n", strerror(errno));
}
size -= did;
buf = (char *)buf - did;
}
}
/* Open a file, and output a tga header for a square image of the given size */
static int open_tga(const char *filename, int size)
{
char c1, c2;
int fd = open(filename, O_CREAT | O_WRONLY | O_TRUNC, S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH);
if (fd == -1) errx(1, "Couldn't open %s for writing\n", filename);
/* Output the header (embedded zeros explicit) */
write(fd, "\0\0\2\0\0\0\0\0\0\0\0\0", 12);
c1 = size%256;
c2 = size/256;
full_write(fd, &c1, 1);
full_write(fd, &c2, 1);
c1 = size%256;
c2 = size/256;
full_write(fd, &c1, 1);
full_write(fd, &c2, 1);
full_write(fd, "\x18\0", 2);
return fd;
}
Where in the above we are careful, and make sure the writes complete correctly. It is easy to just assume that they always will... but it doesn't hurt very much to write a wrapper function that makes 100% sure.
Finally, we can complete our program. This outputs a file called "image.tga" of size 100x100 pixels with a non-spinning black hole (a=0), with a disk viewed at an inclination of 85 degrees.
#define _GNU_SOURCE
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <err.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
#include <errno.h>
#define N 6
static const double a = 0;
static const int size = 100;
static double inclination = 85;
static double r0, theta0;
static double a2;
static double Rhor, Rmstable, Rdisk;
static double L, kappa;
static double inner_orbit(void)
{
double z1 = 1+cbrt(1-a2)*(cbrt(1+a)+cbrt(1-a));
double z2 = sqrt(3*a2+z1*z1);
return 3+z2-sqrt((3-z1)*(3+z1+2*z2));
}
int main(void)
{
int i, j;
int fd = open_tga("image.tga", size);
unsigned char *buffer = calloc(size, 3);
if (!buffer) errx(1, "Out of memory\n");
r0 = 1000.0;
theta0 = (M_PI/180.0) * inclination;
a2 = a*a;
Rhor = 1.0 + sqrt(1.0-a2) + 1e-5;
Rdisk = 20.0;
Rmstable = inner_orbit();
for (j = 0; j < size; j++)
{
for (i = 0; i < size; i++)
{
fire_ray(&buffer[i * 3], i, j);
}
full_write(fd, buffer, size * 3);
printf("%d\n", j);
}
close(fd);
free(buffer);
return 0;
}
Using MPI
The above code doesn't take too long to run, even when the output image is a thousand pixels on a side. However, more realistic simulations require better models of the physics, and perhaps higher-order scattering, where the initial rays can generate many others. The resulting simulations can take weeks or months on a single machine. So it is important to worry about how to make the simulation parallel.
Using MPI, this isn't too difficult. What we do is split the task into two. Most processes can fire rays, and see where they go. However, we need one process to coordinate everything, and output the results. Thus the main() function will turn into:
#include "mpi.h"
int main(int argc, char **argv)
{
int numprocs;
int rank;
/* Initialize MPI */
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if (numprocs == 1) errx(1, "We need more than one rank\n");
if (!rank)
{
master(numprocs - 1);
}
else
{
slave(numprocs - 1, rank - 1);
}
MPI_Finalize();
return 0;
}
Where the "master" function handles the output from the many "slave" processes. This isn't too hard to do if we split the work by scan-line. The master process can just wait for each individual line, and then write them into the write spots in the output file, using the pwrite() function. With MPI, all that is required is a simple call to MPI_Recv for each line. By using MPI_ANY_SOURCE we don't even have to know which scan line will be done next - the MPI library will handle everything. (We can obtain the source from the MPI_SOURCE field of the status structure.
static void full_pwrite(int fd, void *buf, size_t size, size_t offset)
{
while (size)
{
ssize_t did = pwrite(fd, buf, size, offset);
if (did == -1)
{
if (errno == EINTR) continue;
errx(1, "Error writing file: %s\n", strerror(errno));
}
size -= did;
buf = (char *)buf - did;
offset += did;
}
}
/* The master process handles the IO */
static void master(int numslaves)
{
int fd = open_tga("image.tga", size);
int *row = calloc(numslaves, sizeof(int));
off_t offset;
int i;
MPI_Status status;
unsigned char *buf = calloc(size, 3);
if (!buf || !row) errx(1, "Out of memory\n");
/* Initialize the starting rows every slave will be working on */
for (i = 0; i < numslaves; i++) row[i] = i;
/* Wait for size rows of data */
for (i = 0; i < size; i++)
{
/* Wait to get some data */
MPI_Recv(buf, 3 * size, MPI_BYTE, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &status);
/* Calculate offset within output file from the source of the row */
offset = ((off_t) row[status.MPI_SOURCE - 1]) * size * 3 + 18;
/* Output to the correct spot in the file */
full_pwrite(fd, buf, size * 3, offset);
/* The next row it will be working on */
row[status.MPI_SOURCE - 1] += numslaves;
}
close(fd);
free(buf);
free(row);
}
Slightly trickier is the slave code. The most obvious way of doing it is to use MPI_Send after each scan line is complete. The problem with that is that it has to wait until the send is complete before proceeding. If the master is busy handling some other task, we may need to wait quite a while. This is obviously inefficient.
Instead, we can use two buffers and MPI_Isend to make a non-blocking send. It is highly likely that by the time we complete the next buffer that the previous buffer will be handled. This means that the MPI_Wait call will complete immediately without actually having to wait. The result is an almost perfect overlap of communication and computation.
/* Slave compute processes send their results to the master */
static void slave(int numslaves, int rank)
{
MPI_Request rq = MPI_REQUEST_NULL;
int i, j;
unsigned char *buffer1 = calloc(size, 3);
unsigned char *buffer2 = calloc(size, 3);
unsigned char *buf = buffer1;
if (!buffer1 || !buffer2) errx(1, "Out of memory\n");
r0 = 1000.0;
theta0 = (M_PI/180.0) * inclination;
a2 = a*a;
Rhor = 1.0 + sqrt(1.0-a2) + 1e-5;
Rdisk = 20.0;
Rmstable = inner_orbit();
for (j = rank; j < size; j += numslaves)
{
printf("%d\n", j);
for (i = 0; i < size; i++)
{
fire_ray(&buf[i * 3], i, j);
}
/* Wait for the previous send to complete */
MPI_Wait(&rq, MPI_STATUS_IGNORE);
/* Send the new data */
MPI_Isend(buf, 3 * size, MPI_BYTE, 0, 0, MPI_COMM_WORLD, &rq);
/* Flip the buffer to use */
if (buf == buffer1)
{
buf = buffer2;
}
else
{
buf = buffer1;
}
}
/* Wait for the final send to complete */
MPI_Wait(&rq, MPI_STATUS_IGNORE);
/* Free memory */
free(buffer1);
free(buffer2);
}
The above is highly efficient, but we have glossed over something. Note how the master task has a different amount (and type) of work to do than the slaves. With most MPI implementations, this wouldn't be very well designed. The problem is that they will cause the master process to spin inside MPI_Recv wasting cpu resources. What we would like to have happen is the master thread to sleep, and only awaken when needed. Fortunately, Lockless MPI works that way. The master thread will only spin a short time before sleeping.
The fact that processes can go to sleep in Lockless MPI means that load balancing is much simpler than in other MPI implementations. The scan-line hack works well here, but a more advanced code would require the master thread to hand out work units on an as-needed basis. (Otherwise some processes may finish their work well before others, wasting time at the end.) By using a customized control process that load-balances the others, that is easy to accomplish.
Results
We model a disk ranging from 20 Gravitational Radii from the black hole down to the innermost stable orbit. Firstly, we will look at disks around non-rotating "Schwarzschild" black holes. When viewed from 45 degrees, our disk looks like:
The Gravitational Lensing is a little bit subtle in the image, but it is fairly easy to see that the top of the disk has an apparent shape slightly distorted compared to the bottom. However, the biggest feature is the strange ring in the middle, between the innermost stable orbit and the black hole itself. This actually is an image of the bottom of the disk! What happens is that light rays can half-orbit the black hole, and then come from the bottom along those trajectories. Even more difficult to see is a third-order image where light orbits once completely around the hole. In fact there are an infinite series of these images, each smaller than the last.
Now what happens when we view from 85 degrees:
Here, the lensing of the primary image is extreme. The back part of the disk appears to "jump up" due to the photons being bent towards the hole. The second order image below the plane is also becoming quite prominent. (You can also see how left and right sides are swapped there due to the colouration.)
When we move to a spinning black hole, the effects of space-time curvature become more extreme. (We choose a maximally spun-up "Kerr" hole with a=0.998) Viewing again at 45 degrees yields:
Here the most obvious difference is that the disk is much larger. It approaches closer to the event horizon due to the marginally stable orbit decreasing in radius. (And hides the higher order images in the process.) The effects of inertial-frame dragging are also apparent. The innermost part of the image appears to swirl around the black hole. This is caused by the spinning of the black hole imparting a sideways motion onto the photons. (Remember, the disk still has the same shape as it originally did - it just looks different due to the changed lensing.)
Finally, we look at the 85 degree version around our Kerr black hole:
The asymmetry introduced by the spin is highly obvious in this image. We can also see the bottom of the disk peaking out down below. Even though the object near the black hole is a simple 2D ring, the result is quite complex. This shows how difficult it may be to interpret data from such objects. Inverting all the distortions is a complex challenge. However, those distortions themselves can tell us information about the extreme physics involved.
|
Comments
Nobody said...