User: szilagyi
Date: 2013/11/11 05:15 AM
Modified:
/branches/tapir/src/
NullSHRE_WTSphHarm.F90, NullSHRE_modAnaCoord.F90
Log:
work on diagnostic output
File Changes:
Directory: /branches/tapir/src/
===============================
File [modified]: NullSHRE_WTSphHarm.F90
Delta lines: +77 -97
===================================================================
--- branches/tapir/src/NullSHRE_WTSphHarm.F90 2013-11-09 01:48:36 UTC (rev 25)
+++ branches/tapir/src/NullSHRE_WTSphHarm.F90 2013-11-11 11:15:35 UTC (rev 26)
@@ -4,109 +4,68 @@
#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)
+
+subroutine NullSHRE_ConvertAngVector(&
+ ip, Stereo_Vector, Spherical_Vector)
+ ! f^{theta,phi} = f^{q,p} d(theta,phi)/d(q,p)
use NullSHRE_modGFdef
+ use NullSHRE_modAnaCoord
use NullGrid_Vars, only: lsh
- implicit none
+ CCTK_INT, intent(in) :: ip
+ type (gf2d), dimension (2), intent(in) :: Stereo_Vector
+ type (gf2d), dimension (2), intent(inout) :: Spherical_Vector
- 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
+ CCTK_REAL, dimension(lsh(1),lsh(2),4),target :: tmp
+ type (gf2d), dimension(2,2) :: dspherical_dstereo
- ! this will hold the result
- CCTK_REAL, dimension(lsh(1),lsh(2),4) :: dr_Sph
+ 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)
- 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
+ call wt_dspherical_dstereo(ip, dspherical_dstereo)
- ! ! rB = R ( detg g^rr)^(1/4)
+ Spherical_Vector(1)%d = &
+ Stereo_Vector(1)%d * dspherical_dstereo(1,1)%d + &
+ Stereo_Vector(2)%d * dspherical_dstereo(2,1)%d
- ! 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
+ Spherical_Vector(2)%d = &
+ Stereo_Vector(1)%d * dspherical_dstereo(1,2)%d + &
+ Stereo_Vector(2)%d * dspherical_dstereo(2,2)%d
- ! 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
+end subroutine NullSHRE_ConvertAngVector
- ! 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
+subroutine NullSHRE_ConvertAngOneForm(&
+ ip, Stereo_OneForm, Spherical_OneForm)
+ ! f_{theta,phi} = f_{q,p} d(q,p)/d(theta,phi)
+ use NullSHRE_modGFdef
+ use NullSHRE_modAnaCoord
+ use NullGrid_Vars, only: lsh
+ CCTK_INT, intent(in) :: ip
+ type (gf2d), dimension (2), intent(in) :: Stereo_OneForm
+ type (gf2d), dimension (2), intent(inout) :: Spherical_OneForm
+ CCTK_REAL, dimension(lsh(1),lsh(2),4),target :: tmp
+ type (gf2d), dimension(2,2) :: dstereo_dspherical
- !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 )
+ 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 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
+ Spherical_OneForm(1)%d = &
+ Stereo_OneForm(1)%d * dstereo_dspherical(1,1)%d + &
+ Stereo_OneForm(2)%d * dstereo_dspherical(2,1)%d
+ Spherical_OneForm(2)%d = &
+ Stereo_OneForm(1)%d * dstereo_dspherical(1,2)%d + &
+ Stereo_OneForm(2)%d * dstereo_dspherical(2,2)%d
- do a = 2, 4
- write (*,*) 'diff in dr',a,'is: ', maxval(abs(dr_Sph(:,:,a)-dr0(a)%d))
- end do
+end subroutine NullSHRE_ConvertAngOneForm
- ! 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
implicit none
@@ -127,6 +86,7 @@
use NullInterp
use NullDecomp_IO
use NullSHRE_modVars
+ use NullGrid_Vars, only: ip_n, ip_s
implicit none
@@ -142,6 +102,26 @@
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),4),target :: tmp1
+ CCTK_REAL, dimension(null_lsh(1),null_lsh(2),6),target :: tmp2
+ 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
+
+ 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)
+
+ 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)
+
truncate = (IO_TruncateOutputFiles(cctkGH) .ne. 0) .and. first_time
first_time = .FALSE.
@@ -479,7 +459,9 @@
call NullDecomp_WriteCoefsForRealFunc('elldenom_wt', cctkGH, truncate, &
null_lsh, zeta, elld_n%d, elld_s%d, cctk_time)
- ! xb, rb, r_l, r_t
+ ! eta_ab, eta^ab, eta1_ab -- FIXME
+
+ ! 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)
@@ -492,18 +474,16 @@
call NullDecomp_WriteCoefsForRealFunc('r_t_wt', cctkGH, truncate, &
null_lsh, zeta, dr0_n(4)%d, dr0_s(4)%d, cctk_time)
- 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_ConvertAngOneForm(ip_n, dr0_n(2:3), angular_tensor_1_n)
+ call NullSHRE_ConvertAngOneForm(ip_s, dr0_s(2:3), angular_tensor_1_s)
-! 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_theta_wt', cctkGH, truncate, &
+ null_lsh, zeta, angular_tensor_1_n(1)%d, angular_tensor_1_s(1)%d, cctk_time)
+ call NullDecomp_WriteCoefsForRealFunc('r_phi_wt', cctkGH, truncate, &
+ null_lsh, zeta, angular_tensor_1_n(2)%d, angular_tensor_1_s(2)%d, cctk_time)
+
- call NullDecomp_WriteCoefsForRealFunc('r_R_wt', cctkGH, truncate, &
- null_lsh, zeta, temp_n%d, temp_s%d, cctk_time)
+ ! r_A -- FIXME
! J, beta, U, Q, W
File [modified]: NullSHRE_modAnaCoord.F90
Delta lines: +83 -0
===================================================================
--- branches/tapir/src/NullSHRE_modAnaCoord.F90 2013-11-09 01:48:36 UTC (rev 25)
+++ branches/tapir/src/NullSHRE_modAnaCoord.F90 2013-11-11 11:15:35 UTC (rev 26)
@@ -5,6 +5,89 @@
contains
+ subroutine wt_dspherical_dstereo(ip, dspherical_dstereo)
+ use NullSHRE_modGFdef
+ use NullGrid_Vars, only: lsh
+ CCTK_INT, intent(in) :: ip
+ type(gf2d), dimension (2,2), intent(inout) :: dspherical_dstereo
+ type(gf2d), dimension (2,2) :: dstereo_dspherical
+
+ CCTK_REAL, dimension(lsh(1),lsh(2),4),target :: tmp
+ CCTK_REAL, dimension(lsh(1),lsh(2)) :: det
+ 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)
+
+ det = &
+ dstereo_dspherical(1,1)%d * dstereo_dspherical(2,2)%d - &
+ dstereo_dspherical(1,2)%d * dstereo_dspherical(2,1)%d
+
+ call wt_dstereo_dspherical(ip,dstereo_dspherical)
+
+ dspherical_dstereo(1,1)%d = dstereo_dspherical(2,2)%d / det
+ dspherical_dstereo(1,2)%d = - dstereo_dspherical(1,2)%d / det
+ dspherical_dstereo(2,1)%d = dstereo_dspherical(1,1)%d / det
+ dspherical_dstereo(2,2)%d = - dstereo_dspherical(2,1)%d / det
+
+ end subroutine wt_dspherical_dstereo
+
+ subroutine wt_dstereo_dspherical(ip, dstereo_dspherical)
+ use NullSHRE_modGFdef
+ use NullGrid_Vars
+ CCTK_INT, intent(in) :: ip
+ type(gf2d), dimension (2,2), intent(inout) :: dstereo_dspherical
+ CCTK_INT :: pm
+
+ CCTK_REAL, dimension(lsh(1),lsh(2)) :: &
+ CosTheta, SinTheta, CosPhi, SinPhi, x,y,z
+
+ if(ip.eq.ip_n) then
+ pm = +1
+ else
+ pm = -1
+ end if
+
+ ! coordinates on unit sphere
+ x = 2.d0 * qs / pp
+ y = pm * 2.d0 * ps / pp
+ z = pm * (1.d0 - qs ** 2 - ps ** 2) / pp
+
+ ! this will fail if there is a point on the pole
+ if(mod(gsh(1),2) + mod(gsh(2),2) .eq. 2) then
+ call CCTK_WARN(0, "cannot compute wt_dstereo_dspherical on the pole");
+ end if
+
+ CosTheta = z
+ SinPhi = y/sqrt(x*x+y*y)
+ CosPhi = x/sqrt(x*x+y*y)
+ SinTheta = sqrt(1-z*z)
+
+ ! q=(x/r)/(1+pm*(z/r))
+ ! p=pm*(y/r)/(1+pm*(z/r))
+
+ ! x = r cos(phi)sin(theta)
+ ! y = r sin(phi)sin(theta)
+
+ ! q=cos(phi)sin(theta)/(1+pm*cos(theta))
+ ! p=pm*sin(phi)sin(theta)/(1+pm*cos(theta))
+
+ ! theta, phi:
+
+ ! dq/dtheta
+ dstereo_dspherical(1,1)%d = CosPhi/(CosTheta+pm)
+
+ ! dq/dphi
+ dstereo_dspherical(1,2)%d = -SinPhi*SinTheta/(1+pm*CosTheta)
+
+ ! dp/dtheta
+ dstereo_dspherical(2,1)%d = pm*SinPhi/(pm+CosTheta)
+
+ ! dp/dphi
+ dstereo_dspherical(2,2)%d = pm*CosPhi*SinTheta/(1+pm*CosTheta)
+
+ end subroutine wt_dstereo_dspherical
+
subroutine wt_coord (nn1, nn2, qs, ps, pp, ip, cr, wt) ! eq. (17)
use NullSHRE_modGFdef
