Add the capability to rotate initial data, similar to the existing

boost transformation.

Add spherical Kerr-Schild initial data, i.e., Kerr-Schild data in
which the horizon is a coordinate sphere.

Add the capability to smooth Kerr-Schild data with a parabolic term.

Re-indent param.ccl consistently.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@246 e296648e-0e4f-0410-bd07-d597d9acff87

3 files changed, 197 insertions, 1 deletions

diff --git a/src/metrics/Kerr_KerrSchild.F77 b/src/metrics/Kerr_KerrSchild.F77
index f2fb45a..11b6ce1 100644
--- a/src/metrics/Kerr_KerrSchild.F77
+++ b/src/metrics/Kerr_KerrSchild.F77
@@ -38,6 +38,7 @@ c output arguments
 
 c local variables
       CCTK_REAL boostv, eps, m, a
+      integer   power
 
 c local variables
       CCTK_REAL gamma, t0, z0, x0, y0, rho02, r02, r0, costheta0, 
@@ -54,6 +55,7 @@ C     to the J/m definition used in the code here.
 
       boostv = Kerr_KerrSchild__boost_v
       eps    = Kerr_KerrSchild__epsilon
+      power  = Kerr_KerrSchild__power
       m      = Kerr_KerrSchild__mass
       a      = m*Kerr_KerrSchild__spin
 
@@ -79,7 +81,18 @@ C     Spherical auxiliary coordinate r and angle theta in BH rest frame.
 
       r02 = 0.5d0 * (rho02 - a**2) 
      $     + sqrt(0.25d0 * (rho02 - a**2)**2 + a**2 * z0**2)
-      r0 = (r02**2 + eps**4)**0.25d0
+      r0 = sqrt(r02)
+      if (Kerr_KerrSchild__parabolic .eq. 0) then
+C        Use a power law to avoid the singularity
+         r0 = (r0**power + eps**power)**(1.0d0/power)
+      else
+         if (r0 .lt. eps) then
+            r0 = (eps + r0**2 / eps) / 2
+         end if
+      end if
+C     Another idea:
+C        r0 = r0 + eps * exp(-x/eps)
+
       costheta0 = z0 / r0
       
 C     Coefficient H. Note this transforms as a scalar, so does not carry
diff --git a/src/metrics/Kerr_KerrSchild_spherical.F77 b/src/metrics/Kerr_KerrSchild_spherical.F77
new file mode 100644
index 0000000..417b224
--- /dev/null
+++ b/src/metrics/Kerr_KerrSchild_spherical.F77
@@ -0,0 +1,182 @@
+C     Kerr-Schild form of boosted rotating black hole.
+C     Program g_ab = eta_ab + H l_a l_b, g^ab = eta^ab - H l^a l^b.
+C     Here eta_ab is Minkowski in Cartesian coordinates, H is a scalar,
+C     and l is a null vector.
+C     
+C     The coordinates are distorted, such that the event horizon is
+C     a coordinate sphere.
+C
+C     Author: Erik Schnetter <schnetter@cct.lsu.edu>
+C     This formulation was invented by Nils Dorband <dorband@cct.lsu.edu>
+C
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+
+      subroutine Exact__Kerr_KerrSchild_spherical (
+     $     x, y, z, t,
+     $     gdtt, gdtx, gdty, gdtz, 
+     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
+     $     gutt, gutx, guty, gutz, 
+     $     guxx, guyy, guzz, guxy, guyz, guzx,
+     $     psi, Tmunu_flag)
+
+      implicit none
+
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+c input arguments
+      CCTK_REAL x, y, z, t
+
+c output arguments
+      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
+     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
+     $       gutt, gutx, guty, gutz, 
+     $       guxx, guyy, guzz, guxy, guyz, guzx
+      CCTK_REAL psi
+      LOGICAL   Tmunu_flag
+
+c local variables
+      CCTK_REAL boostv, eps, m, a
+
+c local variables
+      CCTK_REAL t1, x1, y1, z1, rho1,
+     $     gamma, t0, z0, x0, y0, rho02, r02, r0, costheta0, 
+     $     lt0, lx0, ly0, lz0, hh, lt, lx, ly, lz
+
+      CCTK_REAL gd(3,3), gdt(3,3), det, jac(3,3)
+
+      integer i, j, k, l
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+C This is a vacuum spacetime with no cosmological constant
+      Tmunu_flag = .false.
+
+C     Get parameters of the exact solution,
+C     and convert from parameter file spin parameter J/m^2
+C     to the J/m definition used in the code here.
+
+      boostv = Kerr_KerrSchild__boost_v
+      eps    = Kerr_KerrSchild__epsilon
+      m      = Kerr_KerrSchild__mass
+      a      = m*Kerr_KerrSchild__spin
+
+C     Distort the coordinates such that the event horizon is a
+C     coordinate sphere
+      rho1 = sqrt (x**2 + y**2 + z**2)
+      rho1 = (rho1**4 + eps**4) ** 0.25d0
+      t1 = t
+      x1 = x - a * y / rho1
+      y1 = y + a * x / rho1
+      z1 = z
+
+C     Boost factor.
+
+      gamma = 1 / sqrt(1 - boostv**2)
+
+C     Lorentz transform t,x,y,z -> t0,x0,y0,z0. 
+C     t0 is never used, but is here for illustration, and we introduce
+C     x0 and y0 also only for clarity.
+C     Note that z0 = 0 means z = vt for the BH.
+
+      t0 = gamma * ((t1 - Kerr_KerrSchild__t) - boostv * (z1 - Kerr_KerrSchild__z))
+      z0 = gamma * ((z1 - Kerr_KerrSchild__z) - boostv * (t1 - Kerr_KerrSchild__t))
+      x0 = x1 - Kerr_KerrSchild__x
+      y0 = y1 - Kerr_KerrSchild__y
+
+C     Spherical auxiliary coordinate r and angle theta in BH rest frame.
+
+      r0 = rho1
+      costheta0 = z0 / r0
+      
+C     Coefficient H.  Note this transforms as a scalar, so does not carry
+C     the suffix 0.
+      hh = m * r0 / (r0**2 + a**2 * costheta0**2)
+      
+C     Components of l_a in rest frame. Note indices down.
+      lt0 = 1.d0
+      lx0 = (r0 * x0 + a * y0) / (r0**2 + a**2)
+      ly0 = (r0 * y0 - a * x0) / (r0**2 + a**2)
+      lz0 = z0 / r0
+
+C     Now boost it to coordinates x, y, z, t.
+C     This is the reverse Lorentz transformation, but applied
+C     to a one-form, so the sign of boostv is the same as the forward
+C     Lorentz transformation applied to the coordinates.
+
+      lt = gamma * (lt0 - boostv * lz0) 
+      lz = gamma * (lz0 - boostv * lt0) 
+      lx = lx0
+      ly = ly0
+
+C     Down metric.  g_ab = flat_ab + H l_a l_b
+
+      gdtt = - 1.d0 + 2.d0 * hh * lt * lt
+      gdtx = 2.d0 * hh * lt * lx
+      gdty = 2.d0 * hh * lt * ly
+      gdtz = 2.d0 * hh * lt * lz
+
+      gd(1,1) = 1.d0 + 2.d0 * hh * lx * lx
+      gd(2,2) = 1.d0 + 2.d0 * hh * ly * ly
+      gd(3,3) = 1.d0 + 2.d0 * hh * lz * lz
+      gd(1,2) = 2.d0 * hh * lx * ly
+      gd(2,3) = 2.d0 * hh * ly * lz
+      gd(3,1) = 2.d0 * hh * lz * lx
+      gd(2,1) = gd(1,2)
+      gd(3,2) = gd(2,3)
+      gd(1,3) = gd(3,1)
+
+C     Transform the tensor basis back
+      jac(1,1) = 1 + (a*x*y) / (rho1**3)
+      jac(1,2) = - a * (x**2+z**2) / (rho1**3)
+      jac(1,3) = a*y*z / (rho1**3)
+      
+      jac(2,1) = a * (y**2+z**2) / (rho1**3)
+      jac(2,2) = 1 - a*x*y / (rho1**3)
+      jac(2,3) = - a*x*z / (rho1**3)
+      
+      jac(3,1) = 0
+      jac(3,2) = 0
+      jac(3,3) = 1
+
+      do i = 1, 3
+         do j = 1, 3
+            gdt(i,j) = 0
+            do k = 1, 3
+               do l = 1, 3
+                  gdt(i,j) = gdt(i,j) + gd(k,l) * jac(k,i) * jac(l,j)
+               end do
+            end do
+         end do
+      end do
+
+      gdxx = gdt(1,1)
+      gdyy = gdt(2,2)
+      gdzz = gdt(3,3)
+      gdxy = gdt(1,2)
+      gdyz = gdt(2,3)
+      gdzx = gdt(3,1)
+
+C     Up metric.  g^ab = flat^ab - H l^a l^b.
+C     Notice that g^ab = g_ab and l^i = l_i and l^0 = - l_0 in flat spacetime.
+      gutt = - 1.d0 - 2.d0 * hh * lt * lt 
+      gutx = 2.d0 * hh * lt * lx
+      guty = 2.d0 * hh * lt * ly
+      gutz = 2.d0 * hh * lt * lz
+
+      det = gdxx*gdyy*gdzz + 2*gdxy*gdzx*gdyz
+     .    - gdxx*gdyz**2 - gdyy*gdzx**2 - gdzz*gdxy**2
+
+      guxx = (gdyy*gdzz - gdyz**2)/det
+      guyy = (gdxx*gdzz - gdzx**2)/det
+      guzz = (gdxx*gdyy - gdxy**2)/det
+
+      guxy = (gdzx*gdyz - gdzz*gdxy)/det
+      guyz = (gdxy*gdzx - gdxx*gdyz)/det
+      guzx = (gdxy*gdyz - gdyy*gdzx)/det
+
+      end
diff --git a/src/metrics/make.code.defn b/src/metrics/make.code.defn
index 12226e2..1e90de0 100644
--- a/src/metrics/make.code.defn
+++ b/src/metrics/make.code.defn
@@ -19,6 +19,7 @@ SRCS = Minkowski.F77			\
        Schwarzschild_Novikov.F77	\
        Kerr_BoyerLindquist.F77		\
        Kerr_KerrSchild.F77		\
+       Kerr_KerrSchild_spherical.F77	\
        Schwarzschild_Lemaitre.F77	\
        multi_BH.F77			\
        Thorne_fakebinary.F77		\