NumericalPolarSchwarzschild.cpp

Demonstration of computing geodesics fully numerically on a discretized four-dimensional grid using polar coordinates. Of course, the computation is slower than using the analytical schwarzschild metric by some orders of magnitude, but the full numerical metric can easily be extended to arbitrary metrices. The goal of the full numerical integration of the Schwarzschild metric is to be able to recover the Photon Orbit, i.e. to get a geodesic making a complete circle around the Schwarzschild black hole. The Photon Orbit is an instable path - it is already hard to compute it even using the analytic solution, so if it can be done in the full numerical regime, we can assume the problem to be solved sufficiently. Ideally, we would also like to be able to trace multiple orbits, but that might be too hard to achieve. The crucial question of interest is: what integrator is required, what step size, what grid size and what interpolator is needed at minimum to get the Photon Orbit? Can we do it in cartesian coordinates as well?

#include <iostream>

#include <spacetime/Schwarzschild.hpp>
#include <spacetime/Geodesic853.hpp>
#include <spacetime/Acceleration.hpp>
#include <vecal/CartesianChart.hpp>
#include <spacetime/NumericalSpacetime.hpp>
#include <spacetime/IntegrateGeodesic.hpp>

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

using namespace VecAl;
using namespace Traum;
using namespace std;

void PolarSchwarzschild4D()
{
        // Use 4D polar coordinates for computation
typedef Polar4Dd::Metric Tensor;
typedef MultiArray<4, Tensor, Tensor*> TensorField;
typedef NumericalSpacetime<TensorField, Polar4Dd> NumSpacetime_t;

        // grid size
enum { N = 10 };

        // define metric tensor field's storage
Tensor TensorStorage[N*N*N*N];
        // define multidimensional field
TensorField G(TensorStorage, MIndex(N,N,N,N) );

        // initialize tensor field by unit metric
Tensor  InitValue;
        InitValue.set(0);
        InitValue(0,0) = InitValue(1,1) = InitValue(2,2) = InitValue(3,3) = 1;
        G.set( InitValue );

        // set up polar metric
        for(int ti=0; ti<N; ti++)
        for(int ri=0; ri<N; ri++)
        for(int hi=0; hi<N; hi++)
        for(int pi=0; pi<N; pi++)
        {
        double  p = pi * 2 * M_PI / N,
                h = hi *     M_PI / N,
                r = ri *     5.0  / N,
                t = ti *     5.0  / N;

        Tensor&g = G[ MIndex(ti, ri, hi, pi) ];

                g(0,0) = 1;//   1 - 2/r;
                g(1,1) = 1;//1/(1 - 2/r);
                g(2,2) = r*r;
                g(3,3) = r*r * sin(p)*sin(p);
        }

        // Define a spacetime with this tensor field
NumSpacetime_t NumSpacetime(G);

        // Define a geodesic integrator on this spacetime
Geodesic853<NumSpacetime_t> myPath(NumSpacetime); 

        // Define initial conditions for the geodesic
Polar4Dd::Point   q; q = 0,3,0,0;
Polar4Dd::Vector  v; v = 1,0,0,1;

ChristoffelField<MultiArray<4, Tensor, Tensor*> > Chris(G, q);

clock_t start, finish;
        start = clock();
        {for(int is = 0; is<10; is++)
        {
                // path parameter increment
        double  ds = 0.001;

                // coordinate acceleration, \ddot q
        Row<4, double> a;

                Chris.getCoordinateAcceleration(a, v);

                // Euler step
                q += v*ds;
                v += a*ds;

        //      if (is %1000 == 0)
                        cout <<" q="<<q<<"\t v="<<v<<"\t a="<< a << endl;
        }}
        finish = clock();
        printf( "4D SW Geodesic: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );
}

void PolarSchwarzschild()
{
        PolarSchwarzschild4D();
}

---