Skip to content
Snippets Groups Projects
Commit f5b22bfe authored by bernuzzi's avatar bernuzzi
Browse files

Added Kerr-Schild metric and extr.curvatuve in spherical polar coordinates for... 

Added Kerr-Schild metric and extr.curvatuve in spherical polar coordinates for any BH spin (https://arxiv.org/abs/1704.00599). To be further tested.
parent a36a52d5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 193 additions and 14 deletions
AHFpy/KerrSchild.py
+ 193
− 14
View file @ f5b22bfe
......@@ -27,19 +27,26 @@ class KerrSchild():
self.x = x
self.y = y
self.z = z
self.grid = None
self.grid = self.CartGrid()
def CartGrid(self):
"""
Cartesian mesh grid
"""
#if not self.grid == None:
# return self.grid
return np.array(np.meshgrid(self.x,self.y,self.z,
indexing="ij"))
def SphereGrid(self):
"""
Spherical Polar mesh grid
"""
x,y,z = self.CartGrid()
rho = np.sqrt(x**2 + y**2)
r = np.sqrt(x**2 + y**2 + z**2)
th = np.arctan2(rho,z)
ph = np.arctan2(y,x)
return r,th,ph
def Rfun(self):
"""
KS Radial function
......@@ -50,7 +57,7 @@ class KerrSchild():
a2 = self.a**2
R2 = 0.5*( rho2 - a2 + np.sqrt((rho2-a2)**2+4*a2*z**2) )
return np.sqrt(R2)
def Vfun(self):
"""
KS Scalar V
......@@ -60,7 +67,7 @@ class KerrSchild():
x,y,z = self.CartGrid()
R = self.Rfun()
return (-self.M * R**3)/(R**4 + self.a**2 * z**2)
def NullCovec(self):
"""
KS Null covector $l_\mu$
......@@ -85,7 +92,7 @@ class KerrSchild():
"""
V = self.Vfun()
return 1.0/np.sqrt(1.0 - 2.0*V)
def Shift_d(self):
"""
KS Shift down index $\beta_i$
......@@ -100,7 +107,7 @@ class KerrSchild():
for i in range(1,4):
beta_d[i] = 2*V*l[i]
return beta_d
def FlatMetric(self):
"""
Minkowski metric (-,+,+,+)
......@@ -147,22 +154,175 @@ class KerrSchild():
for (i,j) in ij_indexes:
K[(i,j)] = fact1 * (R2 * eta[(i,j)] - fact2 * X[i]*X[j])
return K
def Metric4from3(self,alpha,beta_u,gamma):
"""
Build 4-metric $g_{ab}$ from 3-metric $\gamma_{ij}$,
lapse $\alpha$ and shift $\beta^i$
"""
beta_d = {}
beta_d[0] = 0
beta_d[1] = gamma[(1,1)]*beta_u[1] + gamma[(1,2)]*beta_u[2] + gamma[(1,3)]*beta_u[3]
beta_d[2] = gamma[(1,2)]*beta_u[1] + gamma[(2,2)]*beta_u[2] + gamma[(2,3)]*beta_u[3]
beta_d[3] = gamma[(1,3)]*beta_u[1] + gamma[(2,3)]*beta_u[2] + gamma[(3,3)]*beta_u[3]
g = {}
for (i,j) in ij_indexes:
g[(i,j)] = gamma[(i,j)]
g[(0,0)] = - alpha**2
for i in [1,2,3]:
g[(0,0)] += beta_d[i] * beta_u[i]
g[(0,i)] = beta_d[i]
return g
def Output(self,fname):
def MetricSphericalPolar(self):
"""
Metric $g_{ab}$ in spherical polar coordinates
Appendix A of https://arxiv.org/abs/1704.00599
"""
M = self.M
a = self.a
a2 = a**2
r,th,ph = self.SphereGrid()
r2, costh2, sinth2 = r**2, np.cos(th)**2, np.sin(th)**2
rho2 = r2 + a2 * costh2 # Typo in paper!?
# rho = r2 + a2 * costh2
# rho2 = rho**2
tMr_rho2 = 2*M*r/rho2
alpha = 1.0/(np.sqrt(1.0 + tMr_rho2)) # (A.1)
alpha2 = alpha**2
beta_u = {}
beta_u[0] = 0
beta_u[1] = alpha2 * tMr_rho2 # (A.2)
beta_u[2] = 0 # (A.3)
beta_u[3] = 0 # (A.3)
gamma = {}
gamma[(1,1)] = 1.0 + tMr_rho2 # (A.4)
gamma[(1,2)] = 0 # (A.3)
gamma[(1,3)] = -a*gamma[(1,1)]*sinth2 # (A.5)
gamma[(2,2)] = rho2 # (A.3)
gamma[(2,3)] = 0 # (A.7)
gamma[(3,3)] = (r2 + a2 + tMr_rho2*a2*sinth2) * sinth2 # (A.7)
return self.Metric4from3(alpha,beta_u,gamma)
def ExtrinsicCurvatureSphericalPolar(self):
"""
Extrinsic cuvrvature $K_{ij}$ in spherical polar coordinates
Appendix of https://arxiv.org/abs/1704.00599
"""
M = self.M
a = self.a
a2 = a**2
a4 = a**4
r,th,ph = self.SphereGrid()
costh, sinth = np.cos(th), np.sin(th)
cos2th = np.cos(2*th)
r2, costh2, sinth2 = r**2, costh**2, sinth**2
rho2 = r2 + a2 * costh2 # Typo in paper!?
tMr = 2*M*r
tMr_rho2 = tMr/rho2
A = a2*(cos2th + 1.0) + 2*r2 # (A.8)
B = A + 2*tMr # (A.9)
D = np.sqrt(tMr/(a2*costh2 + r2) + 1.0) # (A.10)
K = {}
for (a,b) in ab_indexes:
K[(a,b)] = 0.0
K[(1,1)] = (D*(A+tMr))/(A**2*B)*(4*M*(a2*cos2th+a2-2*r2)) # (A.11)
K[(1,2)] = D/(A*B)*(4*a2*tMr*sinth*costh) # (A.12)
K[(1,3)] = D/A**2*(-2*a*M*sinth2*(a2*cos2th+a2-2*r2)) # (A.13)
K[(2,2)] = D/B*4*M*r2 # (A.14)
K[(2,3)] = D/(A*B)*(-8*a**3*M*r*sinth**3*costh) # (A.15)
K[(3,3)] = D/(A**2*B)*(tMr*sinth2*( a4*(r-M)*np.cos(4*th) + a4*(M+3*r) \
+ 4*a2*r2*(2*r-M) + 4*a2*r*cos2th*(a2+r*(M+2*r)) + 8*r**5 \
) \
) # (A.16-17)
return K
def SphericalPolar2CartesianTensor02(self,g):
"""
Transform tensor $g_{ab}$ from spherical polar to Cartesian components
"""
x,y,z = self.CartGrid()
rho2 = x**2+y**2
rho = np.sqrt(rho2)
r2 = x**2 + y**2 + z**2
r = np.sqrt(r2)
J = {}
J[(0,0)] = 1.0
J[(0,1)] = J[(1,0)] = 0.0
J[(0,2)] = J[(2,0)] = 0.0
J[(0,3)] = J[(3,0)] = 0.0
J[(1,1)] = x/r
J[(1,2)] = y/r
J[(1,3)] = z/r
J[(2,1)] = x*z/(r2*rho)
J[(2,2)] = y*z/(r2*rho)
J[(2,3)] = -rho/r2
J[(3,1)] = -y/(rho2)
J[(3,2)] = x/(rho2)
J[(3,3)] = 0.0
# r,th,ph = self.SphereGrid()
# costh, sinth = np.cos(th), np.sin(th)
# cosph, sinph = np.cos(ph), np.sin(ph)
# J = {}
# J[(0,0)] = 1.0
# J[(0,1)] = J[(1,0)] = 0.0
# J[(0,2)] = J[(2,0)] = 0.0
# J[(0,3)] = J[(3,0)] = 0.0
# J[(1,1)] = sinth*cosph
# J[(1,2)] = r*costh*cosph
# J[(1,3)] = -r*sinth*sinph
# J[(2,1)] = sinth*sinph
# J[(2,2)] = r*costh*sinph
# J[(2,3)] = r*sinth*cosph
# J[(3,1)] = costh
# J[(3,2)] = -r*sinth
# J[(3,3)] = 0.0
# transform to new components
# sort tuple to account for symmetry
gp = {}
for b in range(4):
for a in range(4):
ab_sym = tuple(sorted((a,b)))
gp[ab_sym] = np.zeros_like(r)
for d in range(4):
for c in range(4):
gp[ab_sym] += J[(c,a)] * J[(d,b)] * g[tuple(sorted((c,d)))]
#gp[ab_sym] += J[(a,c)] * J[(b,d)] * g[tuple(sorted((c,d)))]
return gp
def Calculate(self,comp_from='Cartesian'):
"""
Calculate the $g_{ab}$ and $K_{ab}$ in Cartesian components
from either Cartesian or spherical polar expressions
"""
points = self.x, self.y, self.z
if comp_from == 'SphericalPolar':
g = self.MetricSphericalPolar()
K = self.ExtrinsicCurvatureSphericalPolar()
g = self.SphericalPolar2CartesianTensor02(g)
K = self.SphericalPolar2CartesianTensor02(K)
else:
g = self.Metric()
K = self.ExtrinsicCurvature()
return g,K
def Output(self,fname,comp_from='Cartesian'):
"""
Write HDF5 data
"""
g = self.Metric()
K = self.ExtrinsicCurvature()
points = self.x, self.y, self.z
g,K = self.Calculate(comp_from=comp_from)
return write_h5_dset(fname,
time = 0, points = points,
g = g, K = K)
if __name__ == "__main__":
dfile = "kerrschild.h5"
ks = KerrSchild()
ks.Output(dfile)
......@@ -229,3 +389,22 @@ if __name__ == "__main__":
ax.grid()
fig.savefig("kerrschild_alpha.png")
plt.show()
# test data from spherical polar computation
# --------------
dfile = "kerrschild_sp.h5"
ks = KerrSchild()
ks.Output(dfile, comp_from="SphericalPolar")
spdata = read_h5_dset(dfile)
print(spdata.keys())
iz = np.shape(spdata['gxx'])[0]//2
fig, ax = plt.subplots()
ax.contourf(spdata['coord_x'], spdata['coord_y'], spdata['gxy'][:,:,iz])
ax.set(title=r"$g_{xx}$")
ax.grid()
plt.show()
np.testing.assert_allclose(data['gxx'],spdata['gxx'],
rtol=1e-6, atol=0)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment