clean up

src/PointwiseFunctions/GeneralRelativity/CMakeLists.txt

@@ -26,7 +26,6 @@ spectre_target_sources(
   SpacetimeNormalOneForm.cpp
   SpacetimeNormalVector.cpp
   SpatialMetric.cpp
-  SphericalKerrSchildCoords.cpp
   TimeDerivativeOfSpacetimeMetric.cpp
   TimeDerivativeOfSpatialMetric.cpp
   WeylElectric.cpp
@@ -58,7 +57,6 @@ spectre_target_headers(
   SpacetimeNormalOneForm.hpp
   SpacetimeNormalVector.hpp
   SpatialMetric.hpp
-  SphericalKerrSchildCoords.hpp
   Tags.hpp
   TagsDeclarations.hpp
   TimeDerivativeOfSpacetimeMetric.hpp

tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py

@@ -166,7 +166,6 @@ def magnetic_field(
             polytropic_exponent,
         )
     )
-    # FIXME: temporary x in Kerr Schild coordinates.
     x_ks = [0.0, 0.0, 0.0]
 
     if rho > threshold_rho:
@@ -248,7 +247,7 @@ def magnetic_field(
             result, x_ks, bh_mass, bh_dimless_spin
         )
         inv_jac = fm_disk.inverse_jacobian_matrix(x, spin_a)
-        result = np.matmul(inv_jac.T, result)
+        result = inv_jac @ result
     return result
 
 

tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.py

@@ -9,18 +9,18 @@ def delta(r_sqrd, m, a):
 
 
 def transformation_matrix(x, a):
-    # Coordinate transformation matrix for KS to Spherical KS.
+    # Coordinate transformation matrix from SKS to KS.
     # Returns the P^j_ibar from equation 10 of Spherical-KS documentation.
     # Note, this assuemes the black hole spin is along z-direction.
     r_sqrd = boyer_lindquist_r_sqrd(x)
     r = np.sqrt(r_sqrd)
     rho = np.sqrt(r_sqrd + a**2)
-    P_ij = np.diag([rho / r, rho / r, 1.0])
-    return P_ij
+    P_ji = np.diag([rho / r, rho / r, 1.0])
+    return P_ji
 
 
 def jacobian_matrix(x, a):
-    # Jacbian matrix for coordinate transformation from KS to Spherical KS.
+    # Jacbian matrix for coordinate transformation from SKS to KS.
     # Returns T^i_jbar from equation 16 of Spherical-KS documentation.
     # Note, this assuemes the black hole spin is along z-direction.
     P_ij = transformation_matrix(x, a)
@@ -29,10 +29,10 @@ def jacobian_matrix(x, a):
     r = np.sqrt(r_sqrd)
     rho = np.sqrt(r_sqrd + a**2)
 
-    F_ij = (-1.0 / (rho * r**3)) * np.diag([a**2, a**2, 0.0])
+    F_ik = (-1.0 / (rho * r**3)) * np.diag([a**2, a**2, 0.0])
     x_vec = np.array(x).T
 
-    T_ij = P_ij + np.outer(np.matmul(F_ij, x_vec), x_vec.T)
+    T_ij = P_ij + np.outer(np.matmul(F_ik, x_vec), x_vec.T)
 
     return T_ij.T
 
@@ -44,8 +44,7 @@ def inverse_jacobian_matrix(x, a):
     # Note, this assuemes the black hole spin is along z-direction.
     T_ij = jacobian_matrix(x, a)
     S_ij = np.linalg.inv(T_ij)
-
-    return S_ij.T
+    return S_ij
 
 
 def sigma(r_sqrd, sin_theta_sqrd, a):
@@ -91,7 +90,7 @@ def kerr_schild_h(x, m, a):
     return m * r_sqrd * np.sqrt(r_sqrd) / (r_sqrd**2 + a**2 * x[2] ** 2)
 
 
-def kerr_schild_spatial_null_form(x, m, a):
+def kerr_schild_spatial_null_vec(x, m, a):
     r_sqrd = boyer_lindquist_r_sqrd(x)
     r = np.sqrt(r_sqrd)
     rho = np.sqrt(r_sqrd + a**2)
@@ -107,30 +106,20 @@ def kerr_schild_spatial_null_form(x, m, a):
             (r * x[1] * conv_fac - a * x[0] * conv_fac) * denom,
             x[2] / r,
         ]
-    ).T
+    )
+    return l_vec
+
+
+def sph_kerr_schild_spatial_null_form(x, m, a):
     jac = jacobian_matrix(x, a)
-    return (np.matmul(jac, l_vec.T)).T
+    l_vec = kerr_schild_spatial_null_vec(x, m, a)
+    return jac @ l_vec
 
 
-def kerr_schild_spatial_null_vec(x, m, a):
-    r_sqrd = boyer_lindquist_r_sqrd(x)
-    r = np.sqrt(r_sqrd)
-    rho = np.sqrt(r_sqrd + a**2)
-    denom = 1.0 / rho**2
-    # Again, we assume Spherical-KS coordinates
-    # xbar/r = x/rho, ybar/r = y/rho
-    # where rho^2 = r^2+a^2 and xbar and ybar are Spherical-KS
-    # Thus, we need to stick in the converseion factor of rho/r.
-    conv_fac = rho / r
-    l_vec = np.array(
-        [
-            (r * x[0] * conv_fac + a * x[1] * conv_fac) * denom,
-            (r * x[1] * conv_fac - a * x[0] * conv_fac) * denom,
-            x[2] / r,
-        ]
-    ).T
+def sph_kerr_schild_spatial_null_vec(x, m, a):
     inv_jac = inverse_jacobian_matrix(x, a)
-    return (np.matmul(inv_jac, l_vec)).T
+    l_vec = kerr_schild_spatial_null_vec(x, m, a)
+    return inv_jac.T @ l_vec
 
 
 def kerr_schild_lapse(x, m, a):
@@ -141,23 +130,21 @@ def kerr_schild_lapse(x, m, a):
     )
 
 
-def kerr_schild_shift(x, m, a):
+def sph_kerr_schild_shift(x, m, a):
     null_vector_0 = -1.0
     return (
         -2.0
         * kerr_schild_h(x, m, a)
         * null_vector_0
         * kerr_schild_lapse(x, m, a) ** 2
-    ) * kerr_schild_spatial_null_vec(x, m, a)
+    ) * sph_kerr_schild_spatial_null_vec(x, m, a)
 
 
-def kerr_schild_spatial_metric(x, m, a):
+def sph_kerr_schild_spatial_metric(x, m, a):
     prefactor = 2.0 * kerr_schild_h(x, m, a)
-    null_form = kerr_schild_spatial_null_form(x, m, a)
+    null_form = sph_kerr_schild_spatial_null_form(x, m, a)
     T_ij = jacobian_matrix(x, a)
-    result = np.matmul(T_ij, T_ij.T) + prefactor * np.outer(
-        null_form.T, null_form.T
-    )
+    result = T_ij @ T_ij.T + prefactor * np.outer(null_form, null_form)
     return result
 
 
@@ -343,14 +330,14 @@ def spatial_velocity(
                     l, r_sqrd, sin_theta_sqrd, bh_mass, bh_spin_a
                 )
             )
-            transport_velocity_sks = np.matmul(
-                inverse_jacobian_matrix(x, bh_spin_a), transport_velocity_ks.T
+            transport_velocity_sks = (
+                inverse_jacobian_matrix(x, bh_spin_a) @ transport_velocity_ks
             )
             result += (
                 transport_velocity_sks
-                + kerr_schild_shift(x, bh_mass, bh_spin_a)
+                + sph_kerr_schild_shift(x, bh_mass, bh_spin_a)
             ) / kerr_schild_lapse(x, bh_mass, bh_spin_a)
-    return result
+    return np.array(result)
 
 
 def lorentz_factor(
@@ -364,7 +351,7 @@ def lorentz_factor(
     polytropic_exponent,
 ):
     bh_spin_a = bh_mass * bh_dimless_a
-    spatial_metric = kerr_schild_spatial_metric(x, bh_mass, bh_spin_a)
+    spatial_metric = sph_kerr_schild_spatial_metric(x, bh_mass, bh_spin_a)
     velocity = spatial_velocity(
         x,
         t,

tests/Unit/PointwiseFunctions/GeneralRelativity/CMakeLists.txt

@@ -13,7 +13,6 @@ set(LIBRARY_SOURCES
   Test_Psi4.cpp
   Test_Ricci.cpp
   Test_SpacetimeDerivativeOfGothG.cpp
-  Test_SphericalKerrSchildCoords.cpp
   Test_Tags.cpp
   Test_WeylElectric.cpp
   Test_WeylMagnetic.cpp