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

11 Nov 2013

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

    

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

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