[Commits] [svn:einsteintoolkit] NullSHRExtract/branches/tapir/src/ (Rev. 25)
bela at caltech.edu
bela at caltech.edu
Fri Nov 8 19:48:36 CST 2013
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")
More information about the Commits
mailing list