Switching coordinate charts mid-integration with solve_across_charts

I am trying to utilize solve_across_charts to avoid coordinate singularities by switching from spherical to cartesian coordinates when integrating geodesics in the Kerr spacetime. However the code just keeps infinitely running. Following the documentation, my code is (9.3 version),

M = Manifold(4, 'M', latex_name=r'\mathcal{M}', structure='Lorentzian')
X.<t,r,th,ph>=M.chart(r"t r:(2.011,+infinity) th:(0,pi):\theta ph:(0,2*pi):periodic\phi")
C.<T,x,y,z>=M.chart()
X_to_C=X.transition_map(C,[t,r*sin(th)*cos(ph),r*sin(th)*sin(ph), r*cos(th)])
C_to_X=C.transition_map(X,[T,sqrt(x**2+y**2+z**2),atan(sqrt(x**2+y**2)/z),atan(y/x)])
var('m,b,a', domain='real')
m=1
g = M.metric()

rho2 = r^2 + (a*cos(th))^2
Delta = r^2 -2*m*r + a^2
g[0,0] = -(1-2*m*r/rho2)
g[0,3] = -2*a*m*r*sin(th)^2/rho2
g[1,1], g[2,2] = rho2/Delta, rho2
g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/rho2)*sin(th)^2

def g00(m,r,th,a):return g[0,0](m,r,th,a)
def g03(m,r,th,a):return g[0,3](m,r,th,a)
def g11(m,r,th,a):return g[1,1](m,r,th,a)
def g22(m,r,th,a):return g[2,2](m,r,th,a)
def g33(m,r,th,a):return g[3,3](m,r,th,a)
def D(m,r,th,a): return (g03(m,r,th,a))^2-g00(m,r,th,a)*g33(m,r,th,a)

def initial_vector(r0, b,al,ph0=0, E=1, inward=False):
print(b)
t0,th0=0,pi/2
L = -b*E
vth0 = al/(r0*r0)
vt0=(E*g33(m,r0,th0,a)+L*g03(m,r0,th0,a))/(D(m,r0,th0,a))
vr0=sqrt((((1/D(m,r0,th0,a))*((L^2)*g00(m,r0,th0,a)+2*E*L*g03(m,r0,th0,a)+(E^2)*g33(m,r0,th0,a)))
-((vth0)^2)*g22(m,r0,th0,a))/g11(m,r0,th0,a))
if inward:

vr0 = - vr0

vph0 = -(1/D(m,r0,th0,a))*(E*g03(m,r0,th0,a)+L*(g00(m,r0,th0,a)))
p0 = M((t0, r0, th0, ph0), name='p_0')

return M.tangent_space(p0)((vt0, vr0, vth0, vph0), name='v_0')

r0 = 100
th0=pi/2
s = var('s')

ph0=0
v0 = initial_vector(r0, b=0,al=0, ph0=ph0, inward=True)
geod = M.integrated_geodesic(g, (s, 0, 150), v0, across_charts=True)

sol = geod.solve_across_charts(step=0.1, method='odeint', parameters_values={a:0.5},
verbose='False')
interp = geod.interpolate()


I think the issue is on the declaration of the transition maps, however I've tried everything I could think of with no luck. Any help would be much appreciated.

edit retag close merge delete

Sort by » oldest newest most voted

I am not entirely sure how to solve your issue. Indeed If I remove your chart transition, the geodesic declaration (M.integrated_geodesic(...) no longer hangs. My guess is that it is trying to express the metric in the second frame, and fails to invert/simplify something. It does not even start integrating.

But I am wondering what you are trying to achieve with this second chart. None of the systems of coordinates you are using are regular on the horizon. However, Kerr's original coordinates are (see Eq. (3) of this review). Using a coordinate system that is regular on the horizon, it is sufficient to impose $r>0.1$ (or anything small) to avoid the central singularity (no change of variable can save you, since it is a true, physical singularity). solve_across_charts should automatically stop the integration if it exits the chart domain.

If you are interested, you can check this notebook, where I implement basic black-hole ray-tracing step by step. Extending it to the Kerr metric should be straightforward (but very computationally intense, since I make use of the spherical symmetry of Schwarzschild to drastically reduce the computation time).

more

Thanks for your response. My issue is not with the horizon but with the axis of the object. For example in the edge-on case (observer theta=pi/2), a good number of geodesics, with a set of impact parameters (b,a)=(0,a), exit the chart and the integration stops, when it shouldn't.

( 2023-03-09 19:30:40 +0200 )edit

Can you try without solve_across_charts? The integration is supposed to stop when exiting the chart. Are you saying it wrongly detects the edge of the charts? Do you have a minimum working example for the behaviour you describe?

( 2023-03-09 20:18:39 +0200 )edit

For example in the code that is in my post, whatever the impact parameter al maybe, however large, if b=0 meaning that we send a geodesic on the axis of the BH, the integration terminates due to the spherical chart reaching its edge.

( 2023-03-12 19:59:16 +0200 )edit