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; override void Accel(real s, const real *x, const real *v, real *d2x_ds2) { 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 */ } }
1.6.2