NumericalPolarGeodesic2D.cpp

Compute a straight line in polar coordinates (r,phi) in two dimensions to check performance and numerical stability.

Compute a straight line in polar coordinates (r,phi) in two dimensions to check performance and numerical stability.

#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;
 
struct  ExplicitPolarGeodesics2DEquation : IntegratePath853<double>
{
        typedef double real;
 
        void Accel(real s, const real *x, const real *v, real *d2x_ds2) override
        {
                d2x_ds2[0] =   x[0] * v[1] * v[1];
                d2x_ds2[1] = 2/x[0] * v[0] * v[1];
        }
};
 
void ExplicitPolarGeodesics2D()
{
Vector<2,double> q; q = 5,0;
Vector<2,double> v; v = 0,1;
 
ExplicitPolarGeodesics2DEquation PolarPath;
 
        PolarPath.init(2, 0.0, q.ptr(), v.ptr() );
 
        for(int is=0; is<100; is++)
        {
                PolarPath.advance();
                
                if (is %10 == 0)
                {
                double  x = PolarPath.q(0) * sin(PolarPath.q(1)),
                        y = PolarPath.q(0) * cos(PolarPath.q(1));
                        cout << "{" << PolarPath.q(0) << ","<< PolarPath.q(1) << "}" << endl;
                        cout << "\t{x=" << x << ", y=" << y << "}" << endl;
                }
        }
}
 
void PolarGeodesics2D()
{
        //ExplicitPolarGeodesics2D();
 
typedef TangentialSpace<PolarCoordinates<2>, float > Polar2Df;
typedef TangentialSpace<PolarCoordinates<2>, double > Polar2Dd;
typedef Polar2Df myChart;
 
typedef myChart::Metric_t Tensor;
typedef MultiArray<2, Tensor, Tensor*> TensorField;
typedef NumericalSpacetime<TensorField, myChart> NumSpacetime_t;
 
enum { r = PolarCoordinates<2>::r, p = PolarCoordinates<2>::phi }; 
 
        // grid size
enum { N = 10 };
 
        // define metric tensor field's storage
Tensor TensorStorage[N*N];
        // define multidimensional field
TensorField G(TensorStorage, MIndex(N,N) );
 
        // initialize tensor field by unit metric
Tensor  InitValue;
        InitValue.set(0.L);
        InitValue(r,r) = InitValue(p,p) = 1;
        G.set( InitValue );
 
        // set up polar metric
        for(int ri=0; ri<N; ri++)
        for(int pi=0; pi<N; pi++)
        {
        double  P = pi * 2 * M_PI / N,
                R = ri *    5.0  / N;
 
        Tensor&g = G[ MIndex(ri, pi) ];
 
                g(0,0) = 1;
                g(1,1) = R*R;
        }
 
        // Define a spacetime with this tensor field
NumSpacetime_t NumSpacetime(G);
 
        // Define a geodesic integrator on this spacetime
Geodesic853<NumSpacetime_t> myPath(NumSpacetime); 
 
//if (false)
{
        // Define initial conditions for the geodesic
        // Note: q->r = 0 is a coordinate singularity!
myChart::Point   q; 
                  q[q->r] = 5; q[q->phi] = 0;
 
myChart::Vector v; v = 0,1;
        cout <<"Geodesic starts at q="<<q<<"\t v="<<v<< endl;
 
        myPath.restart1(0, q, v);
 
clock_t start, finish;
        start = clock();
        {for(int is = 0; is<10000; is++)
        {
                if (is==0)
                        cout << "Path Accel: " << NumSpacetime(myPath.position(), myPath.velocity() ) << endl;
 
                myPath.advance();
 
        myChart::Point  q = myPath.position();
        myChart::Vector v = myPath.velocity();
 
        double  x = q[q->r] * sin(q[q->phi] ),
                y = q[q->r] * cos(q[q->phi] );
                if (is %1000 == 0)
                {
                        cout <<" q="<<q<<"\tv="<<v;
                        cout << "\t{x=" << x << ", y=" << y << "}" << endl;
                }
        }}
        finish = clock();
        printf( "DOP Geodesic: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );
 
}
{
        // Define initial conditions for the geodesic
        // Note: q->r = 0 is a coordinate singularity!
myChart::Point   q; 
                  q[q->r] = 5; q[q->phi] = 0;
 
myChart::Vector v; v = 0,1;
        cout <<"Geodesic starts at q="<<q<<"\t v="<<v<< endl;
 
typedef ChristoffelField<MultiArray<2, Tensor, Tensor*> > Chris_t;
Chris_t Chris(G, q);
 
clock_t start, finish;
        start = clock();
        {for(int is = 0; is<10000; is++)
        {
                // path parameter increment
        float   ds = 0.001;
 
                // coordinate acceleration, \ddot q
        myChart::Vector a;
 
                Chris.getCoordinateAcceleration(a, v);
 
                if (is==0)
                {
                        cout << "Accel: "<<a<<endl;
                Chris_t::Christoffel_t Gamma;
                        Chris.getLowerChristoffel(Gamma);
                        cout << "Gamma: " << Gamma << endl;
 
                }
 
                // Euler step
                q += v*ds;
                v += a*ds;
 
        double  x = q[q->r] * sin(q[q->phi] ),
                y = q[q->r] * cos(q[q->phi] );
                if (is %1000 == 0)
                {
                        cout <<" q="<<q<<"\tv="<<v;//"\t a="<< a << endl;
                        cout << "\t{x=" << x << ", y=" << y << "}" << endl;
                }
        }}
        finish = clock();
        printf( "2D Geodesic: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );
        /*
        DBL PREC
Geodesic starts at q=(1,0)       v=[0,1]
 q=(1,0.000367496)      v=[-0.000183748,1]      {x=0.000367496, y=1}
 q=(0.967869,0.352351)  v=[-0.16644,0.878467]   {x=0.334017, y=0.908407}
 q=(0.886769,0.630182)  v=[-0.262461,0.627296]  {x=0.522566, y=0.716438}
 q=(0.782168,0.81662)   v=[-0.30019,0.397274]   {x=0.570072, y=0.535541}
 q=(0.669226,0.930705)  v=[-0.312087,0.234651]  {x=0.536746, y=0.399708}
 q=(0.553803,0.996365)  v=[-0.315352,0.130865]  {x=0.464919, y=0.300913}
 q=(0.437733,1.03191)   v=[-0.316143,0.0679784] {x=0.375698, y=0.224635}
 q=(0.321513,1.04955)   v=[-0.316306,0.0314362] {x=0.278817, y=0.1601}
 q=(0.205266,1.05704)   v=[-0.316332,0.0114268] {x=0.178767, y=0.100879}
 q=(0.0890146,1.05924)  v=[-0.316334,0.00200836]        {x=0.0776194, y=0.043575
5}
2D Geodesic: 2.0 seconds
 
        FLT PREC
Geodesic starts at q=(1,0)       v=[0,1]
 q=(1,0.000367496)      v=[-0.000183748,1]      {x=0.000367496, y=1}
 q=(0.967619,0.352353)  v=[-0.168996,0.878504]  {x=0.333932, y=0.908171}
 q=(0.883886,0.630329)  v=[-0.275108,0.628434]  {x=0.520972, y=0.714032}
 q=(0.772552,0.817846)  v=[-0.323736,0.402541]  {x=0.563712, y=0.528267}
 q=(0.649521,0.934973)  v=[-0.342776,0.246054]  {x=0.522594, y=0.385712}
 q=(0.522055,1.00588)   v=[-0.349765,0.147755]  {x=0.440946, y=0.279478}
 q=(0.39298,1.04832)    v=[-0.352269,0.0881573] {x=0.34055, y=0.196109}
 q=(0.263332,1.0736)    v=[-0.353158,0.052478]  {x=0.231449, y=0.125599}
 q=(0.13348,1.08865)    v=[-0.353473,0.0312136] {x=0.118264, y=0.0618924}
 q=(0.00355618,1.0976)  v=[-0.353584,0.0185603] {x=0.00316541, y=0.00162068}
2D Geodesic: 1.5 seconds
        */
}
}