From a29e3c43c4cd68052703cf6f37abdf4c6a5f921a Mon Sep 17 00:00:00 2001 From: schnetter Date: Sat, 21 Apr 2007 02:53:19 +0000 Subject: 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 --- src/boost.F77 | 577 ++++++++++++++++++------------ src/decode_pars.F77 | 2 + src/gauge.F77 | 22 +- src/include/Scalar_CalcTmunu.inc | 2 + src/include/param_defs.inc | 21 +- src/metric.F77 | 9 + src/metrics/Kerr_KerrSchild.F77 | 15 +- src/metrics/Kerr_KerrSchild_spherical.F77 | 182 ++++++++++ src/metrics/make.code.defn | 1 + 9 files changed, 597 insertions(+), 234 deletions(-) create mode 100644 src/metrics/Kerr_KerrSchild_spherical.F77 (limited to 'src') diff --git a/src/boost.F77 b/src/boost.F77 index d81b100..dc26340 100644 --- a/src/boost.F77 +++ b/src/boost.F77 @@ -1,5 +1,7 @@ c This subroutine calculates the 4-metric and its inverse at an event, -c taking into account an optional Lorentz boost. +c taking into account an optional Lorentz boost and an optional rotation. +c The model is first rotated and then boosted, such that the boost is +c applied to the rotated model. c $Header$ c c The coordinates are @@ -9,6 +11,12 @@ c The 4-metrics are c Cgdd(a,b) = Cactus $g_{ab}$ Cguu(a,b) = Cactus $g^{ab}$ c Mgdd(a,b) = Model $g_{ab}$ Mguu(a,b) = Model $g^{ab}$ c +c For a definition of the Euler angles in the conventions used below, see +c http://mathworld.wolfram.com/EulerAngles.html +c Another useful resource may be +c http://en.wikipedia.org/wiki/Euler_angles +c although this uses (on 2006-11-29) different conventions. +c c This file is copyright (c) 2003 by Jonathan Thornburg . c This file is covered by the GNU GPL license; see the files ../README c and ../COPYING for details. @@ -21,7 +29,7 @@ c #include "param_defs.inc" - subroutine Exact__metric( + subroutine Exact__metric( $ decoded_exact_model, $ x, y, z, t, $ gdtt, gdtx, gdty, gdtz, @@ -30,83 +38,94 @@ c $ guxx, guyy, guzz, guxy, guyz, guxz, $ psi) - implicit none - DECLARE_CCTK_FUNCTIONS - DECLARE_CCTK_PARAMETERS + implicit none + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS c input arguments - CCTK_INT decoded_exact_model - CCTK_REAL x, y, z, t + CCTK_INT decoded_exact_model + CCTK_REAL x, y, z, t c output arguments - CCTK_REAL gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz, - $ psi + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ psi c intrinsic functions called - CCTK_REAL sqrt + CCTK_REAL sqrt c static local variables describing Lorentz transformation - logical firstcall - data firstcall /.true./ - CCTK_REAL gamma - CCTK_REAL vv(3), nn(3) - CCTK_REAL parallel(3,3), perp(3,3) - CCTK_REAL Cx_par(3), Cx_perp(3) - CCTK_REAL partial_Mx_wrt_Cx(0:3,0:3) - CCTK_REAL partial_Cx_wrt_Mx(0:3,0:3) - save firstcall - save gamma - save vv, nn - save parallel, perp - save Cx_par, Cx_perp - save partial_Mx_wrt_Cx - save partial_Cx_wrt_Mx + logical firstcall + data firstcall /.true./ + CCTK_REAL gamma + CCTK_REAL vv(3), nn(3) + CCTK_REAL parallel(3,3), perp(3,3) + CCTK_REAL Cx_par(3), Cx_perp(3) + CCTK_REAL partial_Mx_wrt_Cx(0:3,0:3) + CCTK_REAL partial_Cx_wrt_Mx(0:3,0:3) + CCTK_REAL R(0:3,0:3) + save firstcall + save gamma + save vv, nn + save parallel, perp + save Cx_par, Cx_perp + save partial_Mx_wrt_Cx + save partial_Cx_wrt_Mx + save R c coordinates and 4-metric - CCTK_REAL Ct, Cx(3) - CCTK_REAL Cgdd(0:3,0:3), Cguu(0:3,0:3) - CCTK_REAL Mt, Mx(3) - CCTK_REAL Mgdd(0:3,0:3), Mguu(0:3,0:3) + CCTK_REAL Cx(0:3) + CCTK_REAL Cgdd(0:3,0:3), Cguu(0:3,0:3) + CCTK_REAL Mx(0:3) + CCTK_REAL Mgdd(0:3,0:3), Mguu(0:3,0:3) + CCTK_REAL Nx(0:3) + CCTK_REAL Ngdd(0:3,0:3), Nguu(0:3,0:3) c misc temps - CCTK_REAL vnorm, vnormsq - CCTK_REAL delta_ij - CCTK_REAL Cx_par_i, Cx_perp_i - CCTK_REAL vdotCx - CCTK_REAL Cgdd_ab, Cguu_ab - character*100 warn_buffer + CCTK_REAL vnorm, vnormsq + CCTK_REAL delta_ij + CCTK_REAL Cx_par_i, Cx_perp_i + CCTK_REAL vdotCx + CCTK_REAL Cgdd_ab, Cguu_ab + CCTK_REAL cos_phi, sin_phi + CCTK_REAL cos_theta, sin_theta + CCTK_REAL cos_psi, sin_psi + CCTK_REAL R_phi(0:3,0:3), R_theta(0:3,0:3), R_psi(0:3,0:3) + character*1000 warn_buffer c flags, array indices, etc - logical Tmunu_flag - integer i, j - integer Ca, Cb, MA, MB + logical Tmunu_flag + integer i, j, k, l + integer Ca, Cb, MA, MB c constants - integer n - parameter (n = 3) + integer n + parameter (n = 3) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c -c optimized fast-path if no Lorentz boost +c optimized fast-path if no Lorentz boost and no rotation c - if ( (boost_vx .eq. 0.0) - $ .and. (boost_vy .eq. 0.0) - $ .and. (boost_vz .eq. 0.0) ) then - call Exact__metric_for_model( - $ decoded_exact_model, - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz, - $ psi, - $ Tmunu_flag) - return - endif + if ( (boost_vx .eq. 0.0) + $ .and. (boost_vy .eq. 0.0) + $ .and. (boost_vz .eq. 0.0) + $ .and. (rotation_euler_phi .eq. 0.0) + $ .and. (rotation_euler_theta .eq. 0.0) + $ .and. (rotation_euler_psi .eq. 0.0)) then + call Exact__metric_for_model( + $ decoded_exact_model, + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ psi, + $ Tmunu_flag) + return + end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -125,203 +144,321 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c compute Lorentz transformation information on first call c - if (firstcall) then - firstcall = .false. + if (firstcall) then + firstcall = .false. c boost velocity - vv(1) = boost_vx - vv(2) = boost_vy - vv(3) = boost_vz + vv(1) = boost_vx + vv(2) = boost_vy + vv(3) = boost_vz c Lorentz gamma factor, unit vector in direction of boost velocity - vnormsq = 0.0d0 - do 100 i = 1,n - vnormsq = vnormsq + vv(i)*vv(i) -100 continue - gamma = 1.0 / sqrt(1.0 - vnormsq) - vnorm = sqrt(vnormsq) - do 110 i = 1,n - nn(i) = vv(i) / vnorm -110 continue + vnormsq = 0.0 + do i = 1,n + vnormsq = vnormsq + vv(i)*vv(i) + end do + gamma = 1.0 / sqrt(1.0 - vnormsq) + vnorm = sqrt(vnormsq) + if ( (boost_vx .eq. 0.0) + $ .and. (boost_vy .eq. 0.0) + $ .and. (boost_vz .eq. 0.0)) then + nn(1) = 1 + nn(2) = 0 + nn(3) = 0 + else + do i = 1,n + nn(i) = vv(i) / vnorm + end do + end if c projection operators parallel(*,*) and perp(*,*) - do 210 j = 1,n - do 200 i = 1,n - parallel(i,j) = nn(i) * nn(j) - if (i .eq. j) then - delta_ij = 1.0d0 - else - delta_ij = 0.0d0 - endif - perp(i,j) = delta_ij - parallel(i,j) -200 continue -210 continue + do j = 1,n + do i = 1,n + parallel(i,j) = nn(i) * nn(j) + if (i .eq. j) then + delta_ij = 1.0 + else + delta_ij = 0.0 + end if + perp(i,j) = delta_ij - parallel(i,j) + end do + end do c partial derivatives of Model coordinates with respect to Cactus coordinates - partial_Mx_wrt_Cx(0,0) = gamma - do 300 i = 1,n - partial_Mx_wrt_Cx(0,i) = -gamma*vv(i) -300 continue - do 320 i = 1,n - partial_Mx_wrt_Cx(i,0) = -gamma*vv(i) - do 310 j=1,n - partial_Mx_wrt_Cx(i,j) - $ = gamma*parallel(i,j) - $ + perp(i,j) -310 continue -320 continue + partial_Mx_wrt_Cx(0,0) = gamma + do i = 1,n + partial_Mx_wrt_Cx(0,i) = -gamma*vv(i) + end do + do i = 1,n + partial_Mx_wrt_Cx(i,0) = -gamma*vv(i) + do j=1,n + partial_Mx_wrt_Cx(i,j) = gamma*parallel(i,j) + perp(i,j) + end do + end do c partial derivatives of Cactus coordinates with respect to Model coordinates - partial_Cx_wrt_Mx(0,0) = gamma - do 400 i = 1,n - partial_Cx_wrt_Mx(0,i) = + gamma*vv(i) -400 continue - do 420 i = 1,n - partial_Cx_wrt_Mx(i,0) = + gamma*vv(i) - do 410 j=1,n - partial_Cx_wrt_Mx(i,j) - $ = gamma*parallel(i,j) - $ + perp(i,j) -410 continue -420 continue - - endif + partial_Cx_wrt_Mx(0,0) = gamma + do i = 1,n + partial_Cx_wrt_Mx(0,i) = + gamma*vv(i) + end do + do i = 1,n + partial_Cx_wrt_Mx(i,0) = + gamma*vv(i) + do j=1,n + partial_Cx_wrt_Mx(i,j) = gamma*parallel(i,j) + perp(i,j) + end do + end do + +c Sines and cosines of rotation angles + cos_phi = cos (rotation_euler_phi) + sin_phi = sin (rotation_euler_phi) + cos_theta = cos (rotation_euler_theta) + sin_theta = sin (rotation_euler_theta) + cos_psi = cos (rotation_euler_psi) + sin_psi = sin (rotation_euler_psi) + +c Set up individual rotation matrices + R_phi(0,0) = 1 + R_phi(0,1) = 0 + R_phi(0,2) = 0 + R_phi(0,3) = 0 + R_phi(1,0) = 0 + R_phi(1,1) = + cos_phi + R_phi(1,2) = + sin_phi + R_phi(1,3) = 0 + R_phi(2,0) = 0 + R_phi(2,1) = - sin_phi + R_phi(2,2) = + cos_phi + R_phi(2,3) = 0 + R_phi(3,0) = 0 + R_phi(3,1) = 0 + R_phi(3,2) = 0 + R_phi(3,3) = 1 + + R_theta(0,0) = 1 + R_theta(0,1) = 0 + R_theta(0,2) = 0 + R_theta(0,3) = 0 + R_theta(1,0) = 0 + R_theta(1,1) = 1 + R_theta(1,2) = 0 + R_theta(1,3) = 0 + R_theta(2,0) = 0 + R_theta(2,1) = 0 + R_theta(2,2) = + cos_theta + R_theta(2,3) = + sin_theta + R_theta(3,0) = 0 + R_theta(3,1) = 0 + R_theta(3,2) = - sin_theta + R_theta(3,3) = + cos_theta + + R_psi(0,0) = 1 + R_psi(0,1) = 0 + R_psi(0,2) = 0 + R_psi(0,3) = 0 + R_psi(1,0) = 0 + R_psi(1,1) = + cos_psi + R_psi(1,2) = + sin_psi + R_psi(1,3) = 0 + R_psi(2,0) = 0 + R_psi(2,1) = - sin_psi + R_psi(2,2) = + cos_psi + R_psi(2,3) = 0 + R_psi(3,0) = 0 + R_psi(3,1) = 0 + R_psi(3,2) = 0 + R_psi(3,3) = 1 + +c Combine individual rotation matrices + do i = 0,n + do j = 0,n + R(i,j) = 0 + do k = 0,n + do l = 0,n + R(i,j) = R(i,j) + R_psi(i,k) * R_theta(k,l) * R_phi(l,j) + end do + end do + end do + end do + +c Notes that help me (Erik Schnetter) think: +c This considers a rotation with phi=0, theta=pi/2, psi=0. +c Nx(i) = Nx(i) + R(j,i) * Mx(j) +c Nx(1) = - Mx(3) +c Nx(3) = Mx(1) +c Mgxx(1,1) = Ngxx(3,3) +c Mgxx(1,3) = - Ngxx(1,3) +c Mgxx(3,3) = Ngxx(1,1) +c Mgxx(i,j) = R(i,k) R(j,l) Ngxx(k,l) [correct] +c Mgxx(i,j) = R(k,i) R(l,j) Ngxx(k,l) [correct] +c Mbetax(1) = - Nbetax(3) +c Mbetax(3) = Nbetax(1) +c Mbetax(i) + R(i,j) Nbetax(j) [wrong] +c Mbetax(i) + R(j,i) Nbetax(j) [correct] + + end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c compute flat-space components of Cx(*) parallel and perpendicular to vv(*) c - Ct = t - Cx(1) = x - Cx(2) = y - Cx(3) = z - - do 530 i=1,n - Cx_par_i = 0.0d0 - Cx_perp_i = 0.0d0 - do 520 j=1,n - Cx_par_i = Cx_par_i - $ + parallel(i,j)*Cx(j) - Cx_perp_i = Cx_perp_i - $ + perp (i,j)*Cx(j) -520 continue - Cx_par (i) = Cx_par_i - Cx_perp(i) = Cx_perp_i -530 continue + Cx(0) = t + Cx(1) = x + Cx(2) = y + Cx(3) = z + + do i=1,n + Cx_par_i = 0.0 + Cx_perp_i = 0.0 + do j=1,n + Cx_par_i = Cx_par_i + parallel(i,j)*Cx(j) + Cx_perp_i = Cx_perp_i + perp (i,j)*Cx(j) + end do + Cx_par (i) = Cx_par_i + Cx_perp(i) = Cx_perp_i + end do c -c Lorentz-transform the Cactus coordinate to get the Model coordinates +c Lorentz-transform and rotate the Cactus coordinate +c to get the Model coordinates c - vdotCx = 0.0 - do 600 i = 1,n - vdotCx = vdotCx + vv(i)*Cx(i) -600 continue - Mt = gamma * (Ct - vdotCx) - do 610 i=1,n - Mx(i) = gamma * (Cx_par(i) - vv(i)*Ct) - $ + Cx_perp(i) -610 continue +c Boost + vdotCx = 0.0 + do i = 1,n + vdotCx = vdotCx + vv(i)*Cx(i) + end do + + Mx(0) = gamma * (Cx(0) - vdotCx) + do i=1,n + Mx(i) = gamma * (Cx_par(i) - vv(i)*Cx(0)) + Cx_perp(i) + end do + +c Rotation + do i=0,n + Nx(i) = 0 + do j = 0,n + Nx(i) = Nx(i) + R(i,j) * Mx(j) + end do + end do + c c compute the Model 4-metric and inverse 4-metric at the Model coordinates c - call Exact__metric_for_model( - $ decoded_exact_model, - $ Mx(1), Mx(2), Mx(3), Mt, - $ Mgdd(0,0), Mgdd(0,1), Mgdd(0,2), Mgdd(0,3), - $ Mgdd(1,1), Mgdd(2,2), Mgdd(3,3), - $ Mgdd(1,2), Mgdd(2,3), Mgdd(1,3), - $ Mguu(0,0), Mguu(0,1), Mguu(0,2), Mguu(0,3), - $ Mguu(1,1), Mguu(2,2), Mguu(3,3), - $ Mguu(1,2), Mguu(2,3), Mguu(1,3), - $ psi, Tmunu_flag) - - if (Tmunu_flag) then - write (warn_buffer, '(a,i8,a,a)') - $ 'exact_model = ', decoded_exact_model, - $ 'sets the stress-energy tensor', - $ ' ==> we cannot Lorentz-boost it! :(' - call CCTK_WARN(0, warn_buffer) - endif + call Exact__metric_for_model( + $ decoded_exact_model, + $ Nx(1), Nx(2), Nx(3), Nx(0), + $ Ngdd(0,0), Ngdd(0,1), Ngdd(0,2), Ngdd(0,3), + $ Ngdd(1,1), Ngdd(2,2), Ngdd(3,3), + $ Ngdd(1,2), Ngdd(2,3), Ngdd(1,3), + $ Nguu(0,0), Nguu(0,1), Nguu(0,2), Nguu(0,3), + $ Nguu(1,1), Nguu(2,2), Nguu(3,3), + $ Nguu(1,2), Nguu(2,3), Nguu(1,3), + $ psi, Tmunu_flag) + + if (Tmunu_flag) then + write (warn_buffer, '(a,i8,a,a)') + $ 'exact_model = ', decoded_exact_model, + $ 'sets the stress-energy tensor', + $ ' ==> we cannot Lorentz-boost or rotate it!' + call CCTK_WARN(0, warn_buffer) + end if c c symmetrize the Model 4-metric and inverse 4-metric arrays c (the Exact__metric_for_model() call only set the upper triangles) c - Mgdd(1,0) = Mgdd(0,1) - Mgdd(2,0) = Mgdd(0,2) - Mgdd(2,1) = Mgdd(1,2) - Mgdd(3,0) = Mgdd(0,3) - Mgdd(3,1) = Mgdd(1,3) - Mgdd(3,2) = Mgdd(2,3) - - Mguu(1,0) = Mguu(0,1) - Mguu(2,0) = Mguu(0,2) - Mguu(2,1) = Mguu(1,2) - Mguu(3,0) = Mguu(0,3) - Mguu(3,1) = Mguu(1,3) - Mguu(3,2) = Mguu(2,3) + Ngdd(1,0) = Ngdd(0,1) + Ngdd(2,0) = Ngdd(0,2) + Ngdd(2,1) = Ngdd(1,2) + Ngdd(3,0) = Ngdd(0,3) + Ngdd(3,1) = Ngdd(1,3) + Ngdd(3,2) = Ngdd(2,3) + + Nguu(1,0) = Nguu(0,1) + Nguu(2,0) = Nguu(0,2) + Nguu(2,1) = Nguu(1,2) + Nguu(3,0) = Nguu(0,3) + Nguu(3,1) = Nguu(1,3) + Nguu(3,2) = Nguu(2,3) c c tensor-transorm (the upper triangle of) the 4-metric and inverse 4-metric c - do 730 Ca = 0,n - do 720 Cb = Ca,n - Cgdd_ab = 0.0d0 - do 710 Ma = 0,n - do 700 Mb = 0,n - Cgdd_ab = Cgdd_ab - $ + Mgdd(Ma,Mb) - $ * partial_Mx_wrt_Cx(Ma,Ca) - $ * partial_Mx_wrt_Cx(Mb,Cb) -700 continue -710 continue - Cgdd(Ca,Cb) = Cgdd_ab -720 continue -730 continue - - do 830 Ca = 0,n - do 820 Cb = Ca,n - Cguu_ab = 0.0d0 - do 810 Ma = 0,n - do 800 Mb = 0,n - Cguu_ab = Cguu_ab - $ + Mguu(Ma,Mb) - $ * partial_Cx_wrt_Mx(Ca,Ma) - $ * partial_Cx_wrt_Mx(Cb,Mb) -800 continue -810 continue - Cguu(Ca,Cb) = Cguu_ab -820 continue -830 continue + +c Rotations + do i = 0,n + do j = 0,n + Mgdd(i,j) = 0 + Mguu(i,j) = 0 + do k = 0,n + do l = 0,n + Mgdd(i,j) = Mgdd(i,j) + R(k,i) * R(l,j) * Ngdd(k,l) +c The inverse of R is also its transpose. That means that the +c transpose of the inverse, which you would use for g^ij, is just R +c again. + Mguu(i,j) = Mguu(i,j) + R(k,i) * R(l,j) * Nguu(k,l) + end do + end do + end do + end do + +c Boost + do Ca = 0,n + do Cb = Ca,n + Cgdd_ab = 0.0 + do Ma = 0,n + do Mb = 0,n + Cgdd_ab = Cgdd_ab + $ + Mgdd(Ma,Mb) + $ * partial_Mx_wrt_Cx(Ma,Ca) + $ * partial_Mx_wrt_Cx(Mb,Cb) + end do + end do + Cgdd(Ca,Cb) = Cgdd_ab + end do + end do + + do Ca = 0,n + do Cb = Ca,n + Cguu_ab = 0.0 + do Ma = 0,n + do Mb = 0,n + Cguu_ab = Cguu_ab + $ + Mguu(Ma,Mb) + $ * partial_Cx_wrt_Mx(Ca,Ma) + $ * partial_Cx_wrt_Mx(Cb,Mb) + end do + end do + Cguu(Ca,Cb) = Cguu_ab + end do + end do c c unpack the Cactus-coordinates 4-metric and inverse 4-metric c into the corresponding output arguments c - gdtt = Cgdd(0,0) - gdtx = Cgdd(0,1) - gdty = Cgdd(0,2) - gdtz = Cgdd(0,3) - gdxx = Cgdd(1,1) - gdxy = Cgdd(1,2) - gdxz = Cgdd(1,3) - gdyy = Cgdd(2,2) - gdyz = Cgdd(2,3) - gdzz = Cgdd(3,3) - - gutt = Cguu(0,0) - gutx = Cguu(0,1) - guty = Cguu(0,2) - gutz = Cguu(0,3) - guxx = Cguu(1,1) - guxy = Cguu(1,2) - guxz = Cguu(1,3) - guyy = Cguu(2,2) - guyz = Cguu(2,3) - guzz = Cguu(3,3) - - return - end + gdtt = Cgdd(0,0) + gdtx = Cgdd(0,1) + gdty = Cgdd(0,2) + gdtz = Cgdd(0,3) + gdxx = Cgdd(1,1) + gdxy = Cgdd(1,2) + gdxz = Cgdd(1,3) + gdyy = Cgdd(2,2) + gdyz = Cgdd(2,3) + gdzz = Cgdd(3,3) + + gutt = Cguu(0,0) + gutx = Cguu(0,1) + guty = Cguu(0,2) + gutz = Cguu(0,3) + guxx = Cguu(1,1) + guxy = Cguu(1,2) + guxz = Cguu(1,3) + guyy = Cguu(2,2) + guyz = Cguu(2,3) + guzz = Cguu(3,3) + + end diff --git a/src/decode_pars.F77 b/src/decode_pars.F77 index 1b7da46..d51bb4a 100644 --- a/src/decode_pars.F77 +++ b/src/decode_pars.F77 @@ -76,6 +76,8 @@ c black hole spacetimes decoded_exact_model = EXACT__Kerr_BoyerLindquist elseif (CCTK_Equals(exact_model, "Kerr/Kerr-Schild") .ne. 0) then decoded_exact_model = EXACT__Kerr_KerrSchild + elseif (CCTK_Equals(exact_model, "Kerr/Kerr-Schild/spherical") .ne. 0) then + decoded_exact_model = EXACT__Kerr_KerrSchild_spherical elseif (CCTK_Equals(exact_model, "Schwarzschild-Lemaitre") .ne. 0) then decoded_exact_model = EXACT__Schwarzschild_Lemaitre elseif (CCTK_Equals(exact_model, "multi-BH") .ne. 0) then diff --git a/src/gauge.F77 b/src/gauge.F77 index 56088fb..b10bfa1 100644 --- a/src/gauge.F77 +++ b/src/gauge.F77 @@ -19,6 +19,7 @@ C $Header$ integer i,j,k integer nx,ny,nz + CCTK_REAL tt, xx, yy, zz CCTK_REAL gxxjunk, gyyjunk, gzzjunk, $ gxyjunk, gyzjunk, gxzjunk, $ hxxjunk, hyyjunk, hzzjunk, @@ -84,6 +85,11 @@ C do j=1,ny do i=1,nx + tt = cctk_time + xx = x(i,j,k) - cctk_time * shift_add_x + yy = y(i,j,k) - cctk_time * shift_add_y + zz = z(i,j,k) - cctk_time * shift_add_z + C Initialize the psi of exact C (also to tell the models about the conformal_state) if (conformal_state .ne. 0) then @@ -103,7 +109,7 @@ C (also to tell the models about the conformal_state) call Exact__Bona_Masso_data( $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, + $ xx, yy, zz, tt, $ gxxjunk, gyyjunk, gzzjunk, $ gxyjunk, gyzjunk, gxzjunk, $ hxxjunk, hyyjunk, hzzjunk, @@ -147,6 +153,11 @@ C do j=1,ny do i=1,nx + tt = cctk_time + xx = x(i,j,k) - cctk_time * shift_add_x + yy = y(i,j,k) - cctk_time * shift_add_y + zz = z(i,j,k) - cctk_time * shift_add_z + C Initialize the psi of exact C (also to tell the models about the conformal_state) if (conformal_state .ne. 0) then @@ -166,7 +177,7 @@ C (also to tell the models about the conformal_state) call Exact__Bona_Masso_data( $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, + $ xx, yy, zz, tt, $ gxxjunk, gyyjunk, gzzjunk, $ gxyjunk, gyzjunk, gxzjunk, $ hxxjunk, hyyjunk, hzzjunk, @@ -208,6 +219,11 @@ C do j=1,ny do i=1,nx + tt = cctk_time + xx = x(i,j,k) - cctk_time * shift_add_x + yy = y(i,j,k) - cctk_time * shift_add_y + zz = z(i,j,k) - cctk_time * shift_add_z + C Initialize the psi of exact C (also to tell the models about the conformal_state) if (conformal_state .ne. 0) then @@ -227,7 +243,7 @@ C (also to tell the models about the conformal_state) call Exact__Bona_Masso_data( $ decoded_exact_model, - $ x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, + $ xx, yy, zz, tt, $ gxxjunk, gyyjunk, gzzjunk, $ gxyjunk, gyzjunk, gxzjunk, $ hxxjunk, hyyjunk, hzzjunk, diff --git a/src/include/Scalar_CalcTmunu.inc b/src/include/Scalar_CalcTmunu.inc index 346fe09..e4ae8da 100644 --- a/src/include/Scalar_CalcTmunu.inc +++ b/src/include/Scalar_CalcTmunu.inc @@ -61,6 +61,8 @@ c no stress-energy tensor in this model c no stress-energy tensor in this model elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild) then c no stress-energy tensor in this model + elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild_spherical) then +c no stress-energy tensor in this model cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/src/include/param_defs.inc b/src/include/param_defs.inc index 10bbefa..d2374cb 100644 --- a/src/include/param_defs.inc +++ b/src/include/param_defs.inc @@ -38,16 +38,17 @@ c Minkowski spacetime #define EXACT__Minkowski_conf_wave 6 c black hole spacetimes -#define EXACT__Schwarzschild_EF 10 -#define EXACT__Schwarzschild_PG 11 -#define EXACT__Schwarzschild_BL 12 -#define EXACT__Schwarzschild_Novikov 13 -#define EXACT__Kerr_BoyerLindquist 14 -#define EXACT__Kerr_KerrSchild 15 -#define EXACT__Schwarzschild_Lemaitre 16 -#define EXACT__multi_BH 17 -#define EXACT__Alvi 18 -#define EXACT__Thorne_fakebinary 19 +#define EXACT__Schwarzschild_EF 10 +#define EXACT__Schwarzschild_PG 11 +#define EXACT__Schwarzschild_BL 12 +#define EXACT__Schwarzschild_Novikov 13 +#define EXACT__Kerr_BoyerLindquist 14 +#define EXACT__Kerr_KerrSchild 15 +#define EXACT__Kerr_KerrSchild_spherical 16 +#define EXACT__Schwarzschild_Lemaitre 17 +#define EXACT__multi_BH 18 +#define EXACT__Alvi 19 +#define EXACT__Thorne_fakebinary 20 c cosmological spacetimes #define EXACT__Lemaitre 50 diff --git a/src/metric.F77 b/src/metric.F77 index a93ea12..f28c948 100644 --- a/src/metric.F77 +++ b/src/metric.F77 @@ -152,6 +152,15 @@ c $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, $ gutt, gutx, guty, gutz, $ guxx, guyy, guzz, guxy, guyz, guxz, + $ psi, Tmunu_flag) + + elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild_spherical) then + call Exact__Kerr_KerrSchild_spherical( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, $ psi, Tmunu_flag) elseif (decoded_exact_model .eq. EXACT__Schwarzschild_Lemaitre) then 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 +C This formulation was invented by Nils Dorband +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 \ -- cgit v1.2.3