Now the 2d grid array for the generators should work.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@166 2a26948c-0e4f-0410-aee8-f1d3e353619c

2 files changed, 44 insertions, 82 deletions

diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90
index 3dee327..ab19d29 100644
--- a/src/EHFinder_Generator_Sources.F90
+++ b/src/EHFinder_Generator_Sources.F90
@@ -194,18 +194,7 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
         dxg(i,l) = - betaxg(i) + ssign * factor * dfux
         dyg(i,l) = - betayg(i) + ssign * factor * dfuy
         dzg(i,l) = - betazg(i) + ssign * factor * dfuz
-        if ( i == 5 ) then
-          print*,'Debug : ', xg(i,l), yg(i,l), zg(i,l), &
-                             dxg(i,l), dyg(i,l), dzg(i,l)
-            print*,'lapse : ', alp2
-            print*,'shift : ', betaxg(i), betayg(i), betazg(i)
-            print*,'metric : ', gxxg(i), gxyg(i), gxzg(i), &
-                                gyyg(i), gyzg(i), gzzg(i)
-            print*,'imetric : ', guxx, guxy, guxz, guyy, guyz, guzz
-            print*,'df : ', dfxg(i), dfyg(i), dfzg(i)
-            print*,'dfu : ', dfux, dfuy, dfuz
-            print*,'factor : ', factor
-        end if
+
       end do
     else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
       call Util_TableSetIntArray ( status, table_handle, 14, &
@@ -306,8 +295,6 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS)
   CCTK_INT, dimension(14) :: out_types
   CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign
   CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
-  CCTK_REAL, dimension(:,:), allocatable :: tmp, postmp
-  CCTK_INT, dimension(1) :: shape
 
 ! Set the sign correctly depending on the surface direction.
   if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one
@@ -355,12 +342,6 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS)
     call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
   end if
 
-  shape(1) = lsh(1)*lsh(2)
-  allocate ( tmp(lsh(1)*lsh(2),14), postmp(lsh(1)*lsh(2),3) )
-
-  do l = 1, 14
-    out_arrays(l) = CCTK_PointerTo(tmp(:,l))
-  end do
   out_arrays(1) = CCTK_PointerTo(alpg2)
   out_arrays(2) = CCTK_PointerTo(betaxg2)
   out_arrays(3) = CCTK_PointerTo(betayg2)
@@ -376,18 +357,12 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS)
   out_arrays(13) = CCTK_PointerTo(dfzg2)
   out_arrays(14) = CCTK_PointerTo(psig2)
 
-!! Set the pointers to the output arrays.
+! Set the pointers to the output arrays.
 
 ! Loop over the level sets
   do l = 1, eh_number_level_sets
 
 !   Set the pointers to the points to be interpolated to.
-!    postmp(:,1) = reshape(xg2(:,:,l),shape)
-!    postmp(:,2) = reshape(yg2(:,:,l),shape)
-!    postmp(:,3) = reshape(zg2(:,:,l),shape)
-!    interp_coords(1) = CCTK_PointerTo(postmp(:,1))
-!    interp_coords(2) = CCTK_PointerTo(postmp(:,2))
-!    interp_coords(3) = CCTK_PointerTo(postmp(:,3))
 
     interp_coords(1) = CCTK_PointerTo(xg2(1,1,l))
     interp_coords(2) = CCTK_PointerTo(yg2(1,1,l))
@@ -442,20 +417,6 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS)
         call CCTK_INFO ( 'Interpolation failed.' )
       end if
 
-!      alpg2 = reshape(tmp(:,1),lsh)
-!      betaxg2 = reshape(tmp(:,2),lsh)
-!      betayg2 = reshape(tmp(:,3),lsh)
-!      betazg2 = reshape(tmp(:,4),lsh)
-!      gxxg2 = reshape(tmp(:,5),lsh)
-!      gxyg2 = reshape(tmp(:,6),lsh)
-!      gxzg2 = reshape(tmp(:,7),lsh)
-!      gyyg2 = reshape(tmp(:,8),lsh)
-!      gyzg2 = reshape(tmp(:,9),lsh)
-!      gzzg2 = reshape(tmp(:,10),lsh)
-!      dfxg2 = reshape(tmp(:,11),lsh)
-!      dfyg2 = reshape(tmp(:,12),lsh)
-!      dfzg2 = reshape(tmp(:,13),lsh)
-
 !     For each point on this processor calculate the right hand side of the
 !     characteristic evolution equation.
       do j = 1, lsh(2)
@@ -494,19 +455,7 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS)
           dxg2(i,j,l) = - betaxg2(i,j) + ssign * factor * dfux
           dyg2(i,j,l) = - betayg2(i,j) + ssign * factor * dfuy
           dzg2(i,j,l) = - betazg2(i,j) + ssign * factor * dfuz
-          if ( j == 1 .and. i == 5 ) then
-              
-            print*,'debug2 : ', xg2(i,j,l), yg2(i,j,l), zg2(i,j,l), &
-                                dxg2(i,j,l), dyg2(i,j,l), dzg2(i,j,l)
-            print*,'lapse : ', alp2
-            print*,'shift : ', betaxg2(i,j), betayg2(i,j), betazg2(i,j)
-            print*,'metric : ', gxxg2(i,j), gxyg2(i,j), gxzg2(i,j), &
-                                gyyg2(i,j), gyzg2(i,j), gzzg2(i,j)
-            print*,'imetric : ', guxx, guxy, guxz, guyy, guyz, guzz
-            print*,'df : ', dfxg2(i,j), dfyg2(i,j), dfzg2(i,j)
-            print*,'dfu : ', dfux, dfuy, dfuz
-            print*,'factor : ', factor
-          end if
+
         end do
       end do
     else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
@@ -534,21 +483,6 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS)
         call CCTK_INFO ( 'Interpolation failed.' )
       end if
 
-!      alpg2 = reshape(tmp(:,1),lsh)
-!      betaxg2 = reshape(tmp(:,2),lsh)
-!      betayg2 = reshape(tmp(:,3),lsh)
-!      betazg2 = reshape(tmp(:,4),lsh)
-!      gxxg2 = reshape(tmp(:,5),lsh)
-!      gxyg2 = reshape(tmp(:,6),lsh)
-!      gxzg2 = reshape(tmp(:,7),lsh)
-!      gyyg2 = reshape(tmp(:,8),lsh)
-!      gyzg2 = reshape(tmp(:,9),lsh)
-!      gzzg2 = reshape(tmp(:,10),lsh)
-!      dfxg2 = reshape(tmp(:,11),lsh)
-!      dfyg2 = reshape(tmp(:,12),lsh)
-!      dfzg2 = reshape(tmp(:,13),lsh)
-!      psig2 = reshape(tmp(:,14),lsh)
-
       do j = 1, lsh(2)
         do i = 1, lsh(1)
           alp2 = alpg2(i,j)**2
@@ -590,6 +524,5 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS)
     end if
   end do
 
-  deallocate ( tmp, postmp )
   return
 end subroutine EHFinder_Generator_Sources_2D
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index c9dd60c..d9aa13e 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -23,6 +23,7 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
   CCTK_REAL :: last_time
   CCTK_REAL :: theta, dtheta, thetamin, thetamax, r_el
   CCTK_REAL :: phi, dphi, phimin, phimax
+  CCTK_REAL :: costh, sinth, cosph, sinph
   CCTK_INT, dimension(1) :: lsh, lbnd
   CCTK_INT, dimension(2) :: lsh2, lbnd2
 
@@ -115,13 +116,13 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
       if ( number_of_generators_theta .eq. 1 ) then
         theta = half * ( thetamax - thetamin ) + thetamin
       else
-        dtheta = ( thetamax - thetamin ) / ( number_of_generators_theta - 1 )
+        dtheta = ( thetamax - thetamin ) / ( number_of_generators_theta )
       end if
 
       if ( number_of_generators_phi .eq. 1 ) then
         phi = half * ( phimax - phimin ) + phimin
       else
-        dphi = ( phimax - phimin ) / ( number_of_generators_phi - 1 )
+        dphi = ( phimax - phimin ) / ( number_of_generators_phi )
       end if
 
     end if
@@ -140,20 +141,21 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
       if ( evolve_generators .gt. 0 ) then
         if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
           do i = 1, lsh(1)
-            theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
+            theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) + half * dtheta
             xg(i,l) = initial_rad(l) * sin(theta) + translate_x(l)
             yg(i,l) = translate_y(l)
             zg(i,l) = initial_rad(l) * cos(theta) + translate_z(l)
           end do
         else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then
           do j = 1, lsh2(2)
-            phi = phimin + dphi * ( j + lbnd2(2) - 1 )
+            phi = phimin + dphi * ( j + lbnd2(2) - 1 ) + half * dphi
+            cosph = cos(phi); sinph = sin(phi)
             do i = 1, lsh2(1)
               theta = thetamin + dtheta * ( i + lbnd2(1) - 1 )
-              xg2(i,j,l) = initial_rad(l) * sin(theta) * cos(phi) + translate_x(l)
-              yg2(i,j,l) = initial_rad(l) * sin(theta) * sin(phi) + translate_x(l)
-              zg2(i,j,l) = initial_rad(l) * cos(theta) + translate_z(l)
-              print*,xg2(i,j,l),yg2(i,j,l),zg2(i,j,l)
+              costh = cos(theta); sinth = sin(theta)
+              xg2(i,j,l) = initial_rad(l) * sinth * cosph + translate_x(l)
+              yg2(i,j,l) = initial_rad(l) * sinth * sinph + translate_x(l)
+              zg2(i,j,l) = initial_rad(l) * costh + translate_z(l)
             end do
           end do
         end if
@@ -204,16 +206,43 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
         if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
           do i = 1, lsh(1)
             theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
-            r_el = sqrt ( one / ( sin(theta)**2 / initial_a(l)**2 + &
-                           cos(theta)**2 / initial_c(l)**2 ) )
-            xp(1) = r_el * sin(theta) + translate_x(l)
+            costh = cos(theta); sinth = sin(theta)
+            r_el = sqrt ( one / ( sinth**2 / initial_a(l)**2 + &
+                           costh**2 / initial_c(l)**2 ) )
+            xp(1) = r_el * sinth + translate_x(l)
             xp(2) = translate_y(l)
-            xp(3) = r_el * cos(theta) + translate_z(l)
+            xp(3) = r_el * costh + translate_z(l)
             xpt = matmul ( txyz, xp )
             xg(i,l) = xpt(1)
             yg(i,l) = xpt(2)
             zg(i,l) = xpt(3)
           end do
+        else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then
+          do j = 1, lsh2(2)
+            phi = phimin + dphi * ( j + lbnd2(2) - 1 ) + half * dphi
+            cosph = cos(phi); sinph = sin(phi)
+            do i = 1, lsh2(1)
+              theta = thetamin + dtheta * ( i + lbnd2(1) - 1 ) + half * dtheta
+              costh = cos(theta); sinth = sin(theta)
+              r_el = sqrt ( one / &
+                            ( sinth**2*cosph**2 / initial_a(l)**2 + &
+                              sinth**2*sinph**2 / initial_b(l)**2 + &
+                              costh**2 / initial_c(l)**2 ) )
+              xp(1) = r_el * sinth * cosph + translate_x(l)
+              xp(2) = r_el * sinth * sinph + translate_y(l)
+              xp(3) = r_el * costh + translate_z(l)
+              xpt = matmul ( txyz, xp )
+              xg2(i,j,l) = xpt(1)
+              yg2(i,j,l) = xpt(2)
+              zg2(i,j,l) = xpt(3) 
+            end do
+          end do
+!          print*,xg2
+!          print*
+!          print*,yg2
+!          print*
+!          print*,zg2
+!          stop
         end if
       end if
     end if