User: szilagyi Date: 2013/11/11 02:00 PM

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

Log: debugging work

File Changes:

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

File [modified]: NullSHRE_WTSphHarm.F90 Delta lines: +177 -10 =================================================================== --- branches/tapir/src/NullSHRE_WTSphHarm.F90 2013-11-11 11:15:35 UTC (rev 26) +++ branches/tapir/src/NullSHRE_WTSphHarm.F90 2013-11-11 20:00:57 UTC (rev 27) @@ -5,6 +5,84 @@ #include "cctk_Functions.h"

+subroutine NullSHRE_ConvertCoVariantSymTensor(& + ip, Stereo_Tensor, Spherical_Tensor) + ! a=(theta,phi), A=(q,p) + ! T_ab = T_AB dx^A/dx^a dx^B/dx^b + use NullSHRE_modGFdef + use NullSHRE_modAnaCoord + use NullGrid_Vars, only: lsh + CCTK_INT, intent(in) :: ip + type (gf2d), dimension (2,2), intent(in) :: Stereo_Tensor + type (gf2d), dimension (2,2), intent(inout) :: Spherical_Tensor + + CCTK_REAL, dimension(lsh(1),lsh(2),4),target :: tmp + type (gf2d), dimension(2,2) :: dstereo_dspherical + + integer i1,j1,i2,j2 + + dstereo_dspherical(1,1)%d => tmp(:,:,1) + dstereo_dspherical(1,2)%d => tmp(:,:,2) + dstereo_dspherical(2,1)%d => tmp(:,:,3) + dstereo_dspherical(2,2)%d => tmp(:,:,4) + + call wt_dstereo_dspherical(ip, dstereo_dspherical) + + do i1=1,2 + do j1=i1,2 + Spherical_Tensor(i1,j1)%d = 0 + do i2=1,2 + do j2=1,2 + Spherical_Tensor(i1,j1)%d = Spherical_Tensor(i1,j1)%d +& + Stereo_Tensor(i2,j2)%d * dstereo_dspherical(i2,i1)%d& + * dstereo_dspherical(j2,j1)%d + end do + end do + end do + end do + +end subroutine NullSHRE_ConvertCoVariantSymTensor + +subroutine NullSHRE_ConvertContraVariantSymTensor(& + ip, Stereo_Tensor, Spherical_Tensor) + ! a=(theta,phi), A=(q,p) + ! T^ab = T^AB dx^a/dx^A dx^b/dx^B + use NullSHRE_modGFdef + use NullSHRE_modAnaCoord + use NullGrid_Vars, only: lsh + CCTK_INT, intent(in) :: ip + type (gf2d), dimension (2,2), intent(in) :: Stereo_Tensor + type (gf2d), dimension (2,2), intent(inout) :: Spherical_Tensor + + CCTK_REAL, dimension(lsh(1),lsh(2),4),target :: tmp + type (gf2d), dimension(2,2) :: dspherical_dstereo + + integer i1,j1,i2,j2 + + dspherical_dstereo(1,1)%d => tmp(:,:,1) + dspherical_dstereo(1,2)%d => tmp(:,:,2) + dspherical_dstereo(2,1)%d => tmp(:,:,3) + dspherical_dstereo(2,2)%d => tmp(:,:,4) + + call wt_dspherical_dstereo(ip, dspherical_dstereo) + + do i1=1,2 + do j1=i1,2 + Spherical_Tensor(i1,j1)%d = 0 + do i2=1,2 + do j2=1,2 + Spherical_Tensor(i1,j1)%d = Spherical_Tensor(i1,j1)%d +& + Stereo_Tensor(i2,j2)%d * dspherical_dstereo(i2,i1)%d& + * dspherical_dstereo(j2,j1)%d + end do + end do + end do + end do + +end subroutine NullSHRE_ConvertContraVariantSymTensor + + + subroutine NullSHRE_ConvertAngVector(& ip, Stereo_Vector, Spherical_Vector) ! f^{theta,phi} = f^{q,p} d(theta,phi)/d(q,p) @@ -103,25 +181,60 @@ CCTK_REAL, dimension(null_lsh(1), null_lsh(2)) :: fx_s,fy_s,fz_s

CCTK_REAL, dimension(null_lsh(1),null_lsh(2),4),target :: tmp1 - CCTK_REAL, dimension(null_lsh(1),null_lsh(2),6),target :: tmp2 + CCTK_REAL, dimension(null_lsh(1),null_lsh(2),6),target :: & + tmp2, tmp3, tmp4, tmp5 type (gf2d), dimension(2) :: angular_tensor_1_n, angular_tensor_1_s - type (gf2d), dimension(2,2) :: angular_tensor_2s_s, angular_tensor_2s_n + type (gf2d), dimension(2,2) :: eta0_sph_s, eta0_sph_n + type (gf2d), dimension(2,2) :: eta1_sph_s, eta1_sph_n + type (gf2d), dimension(2,2) :: etaup0_sph_s, etaup0_sph_n + type (gf2d), dimension(2,2) :: etaup1_sph_s, etaup1_sph_n

angular_tensor_1_n(1)%d => tmp1(:,:,1) angular_tensor_1_n(2)%d => tmp1(:,:,2) angular_tensor_1_s(1)%d => tmp1(:,:,3) angular_tensor_1_s(2)%d => tmp1(:,:,4)

- angular_tensor_2s_n(1,1)%d => tmp2(:,:,1) - angular_tensor_2s_n(1,2)%d => tmp2(:,:,2) - angular_tensor_2s_n(2,1)%d => tmp2(:,:,2) - angular_tensor_2s_n(2,2)%d => tmp2(:,:,3) + eta0_sph_n(1,1)%d => tmp2(:,:,1) + eta0_sph_n(1,2)%d => tmp2(:,:,2) + eta0_sph_n(2,1)%d => tmp2(:,:,2) + eta0_sph_n(2,2)%d => tmp2(:,:,3)

- angular_tensor_2s_s(1,1)%d => tmp2(:,:,4) - angular_tensor_2s_s(1,2)%d => tmp2(:,:,5) - angular_tensor_2s_s(2,1)%d => tmp2(:,:,5) - angular_tensor_2s_s(2,2)%d => tmp2(:,:,6) + eta0_sph_s(1,1)%d => tmp2(:,:,4) + eta0_sph_s(1,2)%d => tmp2(:,:,5) + eta0_sph_s(2,1)%d => tmp2(:,:,5) + eta0_sph_s(2,2)%d => tmp2(:,:,6)

+ eta1_sph_n(1,1)%d => tmp3(:,:,1) + eta1_sph_n(1,2)%d => tmp3(:,:,2) + eta1_sph_n(2,1)%d => tmp3(:,:,2) + eta1_sph_n(2,2)%d => tmp3(:,:,3) + + eta1_sph_s(1,1)%d => tmp3(:,:,4) + eta1_sph_s(1,2)%d => tmp3(:,:,5) + eta1_sph_s(2,1)%d => tmp3(:,:,5) + eta1_sph_s(2,2)%d => tmp3(:,:,6) + + + etaup0_sph_n(1,1)%d => tmp4(:,:,1) + etaup0_sph_n(1,2)%d => tmp4(:,:,2) + etaup0_sph_n(2,1)%d => tmp4(:,:,2) + etaup0_sph_n(2,2)%d => tmp4(:,:,3) + + etaup0_sph_s(1,1)%d => tmp4(:,:,4) + etaup0_sph_s(1,2)%d => tmp4(:,:,5) + etaup0_sph_s(2,1)%d => tmp4(:,:,5) + etaup0_sph_s(2,2)%d => tmp4(:,:,6) + + etaup1_sph_n(1,1)%d => tmp5(:,:,1) + etaup1_sph_n(1,2)%d => tmp5(:,:,2) + etaup1_sph_n(2,1)%d => tmp5(:,:,2) + etaup1_sph_n(2,2)%d => tmp5(:,:,3) + + etaup1_sph_s(1,1)%d => tmp5(:,:,4) + etaup1_sph_s(1,2)%d => tmp5(:,:,5) + etaup1_sph_s(2,1)%d => tmp5(:,:,5) + etaup1_sph_s(2,2)%d => tmp5(:,:,6) + truncate = (IO_TruncateOutputFiles(cctkGH) .ne. 0) .and. first_time first_time = .FALSE.

@@ -461,6 +574,60 @@

! eta_ab, eta^ab, eta1_ab -- FIXME

+ call NullSHRE_ConvertCoVariantSymTensor(ip_n, eta0_n(2:3,2:3) , eta0_sph_n) + call NullSHRE_ConvertCoVariantSymTensor(ip_s, eta0_s(2:3,2:3) , eta0_sph_s) + call NullSHRE_ConvertCoVariantSymTensor(ip_n, eta1_n(2:3,2:3) , eta1_sph_n) + call NullSHRE_ConvertCoVariantSymTensor(ip_s, eta1_s(2:3,2:3) , eta1_sph_s) + + call NullSHRE_ConvertContraVariantSymTensor(ip_n, etaup0_n(2:3,2:3) , etaup0_sph_n) + call NullSHRE_ConvertContraVariantSymTensor(ip_s, etaup0_s(2:3,2:3) , etaup0_sph_s) + call NullSHRE_ConvertContraVariantSymTensor(ip_n, etaup1_n(2:3,2:3) , etaup1_sph_n) + call NullSHRE_ConvertContraVariantSymTensor(ip_s, etaup1_s(2:3,2:3) , etaup1_sph_s) + + write (*,*) 'stereo identity: 2,2: ', & + maxval(abs( & + -1 + eta0_n(2,2)%d * etaup0_n(2,2)%d + eta0_n(2,3)%d * etaup0_n(3,2)%d )) + + write (*,*) 'stereo identity: 3,3: ', & + maxval(abs( & + -1 + eta0_n(3,2)%d * etaup0_n(2,3)%d + eta0_n(3,3)%d * etaup0_n(3,3)%d )) + + write (*,*) 'stereo identity: 2,3: ', & + maxval(abs( & + eta0_n(2,2)%d * etaup0_n(2,3)%d + eta0_n(2,3)%d * etaup0_n(3,3)%d )) + + write (*,*) 'stereo identity: 3,2: ', & + maxval(abs( & + eta0_n(3,2)%d * etaup0_n(2,2)%d + eta0_n(3,3)%d * etaup0_n(3,2)%d )) + + + + write (*,*) 'spherical identity: 1,1: ', & + maxval(abs( & + -1 + eta0_sph_n(1,1)%d * etaup0_sph_n(1,1)%d + eta0_sph_n(1,2)%d * etaup0_sph_n(2,1)%d )) + + write (*,*) 'spherical identity: 2,2: ', & + maxval(abs( & + -1 + eta0_sph_n(2,1)%d * etaup0_sph_n(1,2)%d + eta0_sph_n(2,2)%d * etaup0_sph_n(2,2)%d )) + + write (*,*) 'spherical identity: 1,2: ', & + maxval(abs( & + eta0_sph_n(1,1)%d * etaup0_sph_n(1,2)%d + eta0_sph_n(1,2)%d * etaup0_sph_n(2,2)%d )) + + write (*,*) 'spherical identity: 2,1: ', & + maxval(abs( & + eta0_sph_n(2,1)%d * etaup0_sph_n(1,1)%d + eta0_sph_n(2,2)%d * etaup0_sph_n(2,1)%d )) + + + + write (*,*) 'rl diff is ', & + maxval(abs( & + -dr0_n(1)%d + 0.25d0*r0_n%d * ( & + + etaup0_sph_n(1,1)%d*eta1_sph_n(1,1)%d& + + etaup0_sph_n(1,2)%d*eta1_sph_n(1,2)%d& + + etaup0_sph_n(2,1)%d*eta1_sph_n(2,1)%d& + + etaup0_sph_n(2,2)%d*eta1_sph_n(2,2)%d ))) + ! xb, rb, r_l, r_t call NullDecomp_WriteCoefsForRealFunc('X_wt', cctkGH, truncate, & null_lsh, zeta, x_wt_n%d, x_wt_s%d, cctk_time)