User: szilagyi Date: 2013/11/08 03:00 PM

Modified: /branches/tapir/src/ NullSHRE_modRadius.F90

Log: work on d(rBondi)/d(rCartesian)

File Changes:

Directory: /branches/tapir/src/ ===============================

File [modified]: NullSHRE_modRadius.F90 Delta lines: +35 -11 =================================================================== --- branches/tapir/src/NullSHRE_modRadius.F90 2013-11-08 19:42:15 UTC (rev 22) +++ branches/tapir/src/NullSHRE_modRadius.F90 2013-11-08 21:00:31 UTC (rev 23) @@ -5,13 +5,13 @@

contains

- subroutine wt_r_long(cr, temp, detg, g3up, j0inv, r0) + subroutine wt_r_long(cr, temp, detg, g3inv, j0inv, r0) use NullSHRE_modGFdef implicit none CCTK_REAL, intent (in) :: cr type (gf2d), intent (in) :: detg type (gf2d), intent (inout) :: temp - type (gf2d), dimension (3,3), intent (in) :: g3up + type (gf2d), dimension (3,3), intent (in) :: g3inv type (gf2d), dimension (4,4), intent (in) :: j0inv type (gf2d), intent (inout) :: r0

@@ -21,7 +21,7 @@ temp%d = 0 do i = 1, 3 do j = 1, 3 - temp%d = temp%d + j0inv(1,i)%d * j0inv(1,j)%d * g3up(i,j)%d + temp%d = temp%d + j0inv(1,i)%d * j0inv(1,j)%d * g3inv(i,j)%d end do end do

@@ -31,28 +31,52 @@

! This is debugging code. given the metric and derivatives, ! will compute drBondi/d(cartesian R) - subroutine wt_dr(g3, detg3, g3inv, dg3_dR, j0inv, r0, dr_dR) + subroutine wt_dr_dR(cr, detg3, g3inv, dg3_dR, j0inv, r0, dr_dR) use NullSHRE_modGFdef + use NullGrid_Vars, only: lsh implicit none

- type (gf2d), dimension (3,3),intent (in) :: g3, g3inv ! g_ij, g^ij + CCTK_REAL, intent (in) :: cr + type (gf2d), dimension (3,3),intent (in) :: g3inv ! g_ij, g^ij type (gf2d), intent (in) :: detg3 ! det(g_ij) type (gf2d), dimension (3,3),intent (in) :: dg3_dR ! g_ij{,r}(cart) + type (gf2d), intent (in) :: r0 ! Bondi r ! x,y,z derivatives of (r, q, p) coordinates : type (gf2d), dimension (3,3),intent (in) :: j0inv - type (gf2d), intent (in) :: r0 + ! this will hold the result type (gf2d), intent (inout) :: dr_dR

+ CCTK_REAL, dimension(lsh(1),lsh(2)) :: detg_R, grr, grr_R, g3invij_R + integer :: i, j, k, l +! ! radial deriv of j0inv: +! ! dj0inv(a,:)/dr = 0 for a=1 +! ! dj0inv(a,:)/dr = -j0inv(a,:)/cr for a=2,3

- ! radial deriv of j0inv: - ! dj0inv(a,:)/dr = 0 for a=1 - ! dj0inv(a,:)/dr = -j0inv(a,:)/cr for a=2,3 +! ! rB = R ( detg g^rr)^(1/4)

- ! rB = R ( detg g^rr)^(1/4) + detg_R = 0 + grr = 0 + grr_R = 0 + g3invij_R = 0

+ do i = 1, 3 + do j = 1, 3 + detg_R = detg_R + g3inv(i,j)%d * dg3_dR(i,j)%d + grr = grr + j0inv(1,i)%d * j0inv(1,j)%d * g3inv(i,j)%d + do k = 1, 3 + do l = 1, 3 + g3invij_R = g3invij_R-g3inv(i,k)%d*g3inv(j,l)%d*dg3_dR(k,l)%d + end do + end do + ! note that partial_R j0inv(1,:) is zero + grr_R = grr_R + j0inv(1,i)%d * j0inv(1,j)%d * g3invij_R + end do + end do

+! dr_dR%d = r0%d / R + 0.25 * R^4 * ( detg g^rr )^(-3/4) / R^3 * ()_,r + dr_dR%d = r0%d/cr + 0.25*cr * (cr/r0%d)**3 * ( detg_R*grr + detg3%d*grr_R )

- end subroutine wt_dr + end subroutine wt_dr_dR

subroutine wt_r(nn1, nn2, pp, eta0, r0) ! the bondi radius eq. (53)