User: szilagyi Date: 2013/11/08 07:48 PM

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

Log: work on alternatives of dr

File Changes:

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

File [modified]: NullSHRE_WTSphHarm.F90 Delta lines: +117 -5 =================================================================== --- branches/tapir/src/NullSHRE_WTSphHarm.F90 2013-11-08 21:55:39 UTC (rev 24) +++ branches/tapir/src/NullSHRE_WTSphHarm.F90 2013-11-09 01:48:36 UTC (rev 25) @@ -4,9 +4,111 @@ #include "cctk_Parameters.h" #include "cctk_Functions.h"

+! This is debugging code. given the metric and derivatives, +! will compute drBondi/d(cartesian R) +subroutine NullSHRE_Calc_dr_Sph( & + cr, detg3, g3inv, dg_Sph, j0inv, r0, dr_dR, dr0, ell) + use NullSHRE_modGFdef + use NullGrid_Vars, only: lsh + implicit none + + 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 (4,4,4),intent (in) :: dg_Sph ! g_ij{,r}(cart) + type (gf2d), intent (in) :: r0 ! Bondi r + ! x,y,z derivatives of (r, q, p) coordinates : + type (gf2d), dimension (4,4), intent (in) :: j0inv + ! this is ell^a and dr/d(lambda,q,p,u) -- used for consistency tests + type (gf2d), dimension (4), intent (in) :: dr0, ell + type (gf2d), intent (inout) :: dr_dR + + ! this will hold the result + CCTK_REAL, dimension(lsh(1),lsh(2),4) :: dr_Sph + + CCTK_REAL, dimension(lsh(1),lsh(2),4) :: ddetg_Sph, dgrr_Sph + CCTK_REAL, dimension(lsh(1),lsh(2)) :: grr, r1 + integer :: i, j, k, l, a + + ! ! rB = R ( detg g^rr)^(1/4) + + ! derivatives of detg + ! d(detg) = detg * ( g^ij d(g_ij) ) + do a = 1, 4 + ddetg_Sph(:,:,a) = 0 + do i = 1, 3 + do j = 1, 3 + ddetg_Sph(:,:,a) = ddetg_Sph(:,:,a)+g3inv(i,j)%d*dg_Sph(i,j,a)%d + end do + end do + ddetg_Sph(:,:,a) = ddetg_Sph(:,:,a) * detg3%d + end do + + ! g^rr = g^ij invJ^r_i invJ^r_j + grr = 0 + do i = 1, 3 + do j = 1, 3 + grr = grr + j0inv(1,i)%d * j0inv(1,j)%d * g3inv(i,j)%d + end do + end do + + ! derivatives of g^rr + ! NOTE: + ! radial deriv of j0inv: + ! dj0inv(a,:)/dr = 0 for a=1 or 4 + ! dj0inv(a,:)/dr = -j0inv(a,:)/cr for a=2,3 + + do a = 1, 4 + dgrr_Sph(:,:,a) = 0 + do i = 1, 3 + do j = 1, 3 + do k = 1, 3 + do l = 1, 3 + dgrr_Sph(:,:,a) = & + dgrr_Sph(:,:,a) - j0inv(1,i)%d * j0inv(1,j)%d * & + g3inv(i,k)%d * g3inv(i,l)%d * dg_Sph(k,l,a)%d + end do + end do + if(a.eq.2 .or. a.eq.3) then + dgrr_Sph(:,:,a) = dgrr_Sph(:,:,a) - 2.0d0 * grr / cr + end if + end do + end do + end do + + + !r0%d = cr * (detg*grr)**0.25d0 + ! dr_dR%d = r0%d/cr + 0.25*cr * (cr/r0%d)**3 * ( detg_R*grr + detg3%d*grr_R ) + + + do a = 1, 4 + dr_Sph(:,:,a) = cr * 0.25d0 * (detg3%d*grr)**(0.25d0-1.0d0) * & + ( ddetg_Sph(:,:,a) * grr + detg3%d*dgrr_Sph(:,:,a)) + if(a.eq.1) dr_Sph(:,:,a) = dr_Sph(:,:,a) + r0%d / cr + end do + + + do a = 2, 4 + write (*,*) 'diff in dr',a,'is: ', maxval(abs(dr_Sph(:,:,a)-dr0(a)%d)) + end do + + ! 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_Sph(:,:,1) + 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)) + + dr_dR%d = dr_Sph(:,:,1) +end subroutine NullSHRE_Calc_dr_Sph + 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 @@ -37,8 +139,8 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS

- 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 + 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

truncate = (IO_TruncateOutputFiles(cctkGH) .ne. 0) .and. first_time first_time = .FALSE. @@ -390,9 +492,19 @@ 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) + call NullSHRE_Calc_dr_Sph( & + cr, detg_n, g3up_n, dg_n, j0inv_n, r0_n, temp_n, dr0_n, ell_n) + call NullSHRE_Calc_dr_Sph( & + cr, detg_s, g3up_s, dg_s, j0inv_s, r0_s, temp_s, dr0_s, ell_s) +! call NullSHRE_Calc_r_R(cr, detg_n, g3up_n, & +! dg_n(1:3,1:3,1), j0inv_n(1:3,1:3), r0_n, temp_n, dr0_n, ell_n)

+! call NullSHRE_Calc_r_R(cr, detg_s, g3up_s, & +! dg_s(1:3,1:3,1), j0inv_s(1:3,1:3), r0_s, temp_s, dr0_s, ell_s) + + call NullDecomp_WriteCoefsForRealFunc('r_R_wt', cctkGH, truncate, & + null_lsh, zeta, temp_n%d, temp_s%d, cctk_time) + ! J, beta, U, Q, W

call NullDecomp_WriteCoefsForComplexFunc('J_wt', cctkGH, truncate, &

File [modified]: NullSHRE_modRadius.F90 Delta lines: +5 -67 =================================================================== --- branches/tapir/src/NullSHRE_modRadius.F90 2013-11-08 21:55:39 UTC (rev 24) +++ branches/tapir/src/NullSHRE_modRadius.F90 2013-11-09 01:48:36 UTC (rev 25) @@ -29,69 +29,6 @@

end subroutine wt_r_long

- ! 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, dr0, ell) - use NullSHRE_modGFdef - use NullGrid_Vars, only: lsh - implicit none - - 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 - ! 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, 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 - -! ! 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 ) - - ! 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) ! the bondi radius eq. (53) use NullSHRE_modGFdef @@ -121,6 +58,7 @@ CCTK_INT :: a, b

character(200) message + CCTK_REAL min_rl

! lambda derivative of bondi r

@@ -132,10 +70,10 @@ end do dr0(1)%d = (r0%d / 4.d0) * temp%d

- if (halt_on_negative_rl.ne.0) then - if (any(dr0(1)%d.lt.0)) call CCTK_WARN(0, "negative value for r1_wt") - else - if (any(dr0(1)%d.lt.rl_min)) then + if (any(dr0(1)%d.lt.rl_min)) then + if (halt_on_negative_rl.ne.0) then + if (any(dr0(1)%d.lt.0)) call CCTK_WARN(0, "negative value for r1_wt") + else write (message,'("r_{,lambda} is less than minimum accepted value: ",g15.5," < ",g15.5)') minval(dr0(1)%d), rl_min call CCTK_INFO(message) call CCTK_INFO("adding artificial correction term to r1_wt")