[Commits] [svn:einsteintoolkit] NullSHRExtract/branches/tapir/src/ (Rev. 23)

[bela at caltech.edu]

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)