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

[bela at caltech.edu]

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

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

Log:
 work on dr(Bondi)/dr(Cartesian)

File Changes:

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

File [modified]: NullSHRE_WTSphHarm.F90
Delta lines: +6 -2
===================================================================
--- branches/tapir/src/NullSHRE_WTSphHarm.F90	2013-11-08 21:00:31 UTC (rev 23)
+++ branches/tapir/src/NullSHRE_WTSphHarm.F90	2013-11-08 21:55:39 UTC (rev 24)
@@ -6,6 +6,7 @@
 
 subroutine NullSHRE_ConvertDerivs(null_lsh, fr, fq, fp, fx, fy, fz, j0inv)
   use NullSHRE_modGFdef
+  use NullSHRE_modRadius, only: wt_dr_dR
   implicit none
   CCTK_INT,  dimension(2),                        intent(in)  :: null_lsh
   CCTK_REAL, dimension(null_lsh(1), null_lsh(2)), intent(in)  :: fr, fq, fp
@@ -36,8 +37,8 @@
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
-  CCTK_REAL, dimension(null_lsh(1), null_lsh(2)) :: fx_n,fy_n,fz_n
-  CCTK_REAL, dimension(null_lsh(1), null_lsh(2)) :: fx_s,fy_s,fz_s
+  CCTK_REAL, dimension(null_lsh(1), null_lsh(2)) :: fx_n,fy_n,fz_n, rR_n
+  CCTK_REAL, dimension(null_lsh(1), null_lsh(2)) :: fx_s,fy_s,fz_s, rR_s
 
   truncate = (IO_TruncateOutputFiles(cctkGH) .ne. 0) .and. first_time
   first_time = .FALSE.
@@ -389,6 +390,9 @@
   call NullDecomp_WriteCoefsForRealFunc('r_t_wt', cctkGH, truncate, &
      null_lsh, zeta, dr0_n(4)%d, dr0_s(4)%d, cctk_time)
 
+   call wt_dr_dR(cr, detg_n, g3up_n, &
+     dg_n(1:3,1:3,1), j0inv_n(1:3,1:3), r0_n, rR_n, dr0_n, ell_n)
+
    ! J, beta, U, Q, W
 
   call NullDecomp_WriteCoefsForComplexFunc('J_wt', cctkGH, truncate, &

File [modified]: NullSHRE_modRadius.F90
Delta lines: +17 -3
===================================================================
--- branches/tapir/src/NullSHRE_modRadius.F90	2013-11-08 21:00:31 UTC (rev 23)
+++ branches/tapir/src/NullSHRE_modRadius.F90	2013-11-08 21:55:39 UTC (rev 24)
@@ -31,7 +31,7 @@
 
   ! This is debugging code.  given the metric and derivatives,
   ! will compute drBondi/d(cartesian R)
-  subroutine wt_dr_dR(cr, detg3, g3inv, dg3_dR, j0inv, r0, dr_dR)
+  subroutine wt_dr_dR(cr, detg3, g3inv, dg3_dR, j0inv, r0, dr_dR, dr0, ell)
     use NullSHRE_modGFdef
     use NullGrid_Vars, only: lsh
     implicit none
@@ -45,9 +45,11 @@
     type (gf2d), dimension (3,3),intent (in)    :: j0inv
     ! this will hold the result
     type (gf2d),                 intent (inout) :: dr_dR
+    ! this is ell^a and dr/d(lambda,q,p,u) -- used for consistency tests
+    type (gf2d), dimension (4),  intent (in)    :: dr0, ell
 
-    CCTK_REAL, dimension(lsh(1),lsh(2)) :: detg_R, grr, grr_R, g3invij_R
-    integer :: i, j, k, l
+    CCTK_REAL, dimension(lsh(1),lsh(2)) :: detg_R, grr, grr_R, g3invij_R, r1
+    integer :: i, j, k, l, a
 !   ! radial deriv of j0inv:
 !   !   dj0inv(a,:)/dr = 0 for a=1
 !   !   dj0inv(a,:)/dr = -j0inv(a,:)/cr for a=2,3
@@ -76,6 +78,18 @@
 !   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 )
 
+    ! build test:
+    ! r_{,lambda} = ell^t r_{,u} + ell^i InvJac^k_i d(rBondi)/dy^k
+
+    r1 = ell(4)%d * dr0(4)%d
+    do i = 1, 3
+       r1 = r1 + ell(i)%d * j0inv(1,i)%d * dr_dR%d
+       do a = 2, 3
+          r1 = r1 + ell(i)%d * j0inv(a,i)%d * dr0(a)%d
+       end do
+    end do
+    write (*,*) 'diff in r1 is ', maxval(abs(r1 - dr0(1)%d))
+
   end subroutine wt_dr_dR
 
   subroutine wt_r(nn1, nn2, pp, eta0, r0)