Added support for finding out how many surfaces are present and then do the

integrals over all surfaces. Still only supports full mode. Only weakly tested.
There is a problem with the outputs in multiprocessor mode. The results of
the integration are stored in 1D grid arrays. However if I use DISTRIB=DEFAULT
(taking care only to store results on the correct processor), the
output fails if the number of processors are bigger than the number of
elements in the grid array (essentially some processors contains a chunk
of length zero). If I on the other hand uses DISTRIB=CONSTANT, it works
but on multiple processors i get multiple copies of the data in the output.


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

2 files changed, 149 insertions, 34 deletions

diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index fe41fd1..41bbac4 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -15,11 +15,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: i, j, k, im, jm
+  CCTK_INT :: i, j, k, im, jm, sn
   CCTK_REAL, parameter :: eps = 1.0d-10
   CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
   CCTK_INT :: interp_handle2
 
+  character(len=30) :: info_message
   CCTK_INT, dimension(4) :: bbox
   CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost, maxl
 
@@ -42,7 +43,10 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 ! the 3D grid functions.
 #include "include/physical_part.h"
 
-  call CCTK_INFO ( 'Finding points on surface' )
+  sn = surface_counter - integrate_counter + 1
+
+  write(info_message,'(a26,i3)') 'Finding points on surface ',sn
+  call CCTK_INFO ( info_message )
 
   ! Find information about the local and global properties of the 2D
   ! Grid arrays.
@@ -77,26 +81,33 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   end if
   
   ! If the user defined center point should be used...
-  if ( use_user_center ) then
+  if ( use_user_center .gt. 0 ) then
     center(1) = center_x
     center(2) = center_y
     center(3) = center_z
   else
 
     ! IF not find an approximation for it by averaging the coodinates that
-    ! are within one grid cell distance from the surface.
+    ! are within one grid cell distance from the surface in question.
     nc_loc = 0; center_loc = zero
     delta = minval ( cctk_delta_space )
     do k = kzl, kzr
       do j = jyl, jyr
         do i = ixl, ixr
           if ( eh_mask(i,j,k) .eq. 0 ) then
-            if ( ( f(i,j,k) * f(i-1,j,k) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i+1,j,k) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i,j-1,k) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i,j+1,k) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i,j,k-1) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i,j,k+1) .lt. zero ) ) then
+            if ( ( ( f(i,j,k) * f(i-1,j,k) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i+1,j,k) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i,j-1,k) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i,j+1,k) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i,j,k-1) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i,j,k+1) .lt. zero ) ) .and. &
+                 ( ( sc(i,j,k) .eq. sn ) .or. &
+                   ( sc(i-1,j,k) .eq. sn ) .or. & 
+                   ( sc(i+1,j,k) .eq. sn ) .or. &
+                   ( sc(i,j-1,k) .eq. sn ) .or. &
+                   ( sc(i,j+1,k) .eq. sn ) .or. &
+                   ( sc(i,j,k-1) .eq. sn ) .or. &
+                   ( sc(i,j,k+1) .eq. sn ) ) ) then
               nc_loc = nc_loc + 1
               center_loc(1) = center_loc(1) + x(i,j,k)
               center_loc(2) = center_loc(2) + y(i,j,k)
@@ -115,7 +126,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                         CCTK_VARIABLE_REAL )
     ncinv = one / nc
     center = ncinv * center
-!    print*,center
+!    print*,'EHFinder_FindSurfaces : ', sn, center
 
     ! For now symmetries are not handled properly
     if ( symx ) center(1) = zero
@@ -183,8 +194,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   ! Get the 3D coordinate system handle.
   call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
   if ( coord_system_handle .lt. 0) then
-    call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. &
- &  Forgot to activate an implementation providing coordinates ??" )
+    call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system.  Forgot to activate an implementation providing coordinates ??" )
   endif
 
   ! Get pointers to the grid arrays containing the interpolation points.
@@ -427,6 +437,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 
   call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" )
 
+  integrate_counter = integrate_counter - 1
 !  print*,'Calculated derivatives'
 !  print*,rsurf
   call Util_TableDestroy ( status, table_handle )
@@ -634,4 +645,7 @@ subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS)
   end if
   call CCTK_INFO(info_message)
 
+  integrate_counter = surface_counter
+
+!  print*,'EHFinder_InfoSurfaces : ',integrate_counter
 end subroutine EHFinder_InfoSurfaces
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
index f82cd54..e463a43 100644
--- a/src/EHFinder_Integrate2.F90
+++ b/src/EHFinder_Integrate2.F90
@@ -72,8 +72,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
   
   call CCTK_InterpHandle ( interp_handle, "Lagrange polynomial interpolation" )
   if ( interp_handle .lt. 0 ) then
-    call CCTK_WARN( 0, "Cannot get handle for interpolation. &
- &  Forgot to activate an implementation providing interpolation operators??" )
+    call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
   end if
 
   call Util_TableCreateFromString ( table_handle, "order=3" )
@@ -83,8 +82,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
   
   call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
   if ( coord_system_handle .lt. 0) then
-    call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. &
- &  Forgot to activate an implementation providing coordinates ??" )
+    call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. Forgot to activate an implementation providing coordinates ??" )
   endif
 
   interp_coords(1) = CCTK_PointerTo(interp_x)
@@ -201,10 +199,14 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: int_var
+  CCTK_INT :: int_var, sn
+  CCTK_REAL :: eh_area_tmp
   character(len=30) :: info_message
+!  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
 
   call CCTK_INFO ( 'Integrating area' )
+
+  sn = surface_counter - integrate_counter
   int_tmp = weights * da
 
   call CCTK_VarIndex ( int_var, "ehfinder::int_tmp" )
@@ -212,10 +214,33 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
     Call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' )
   end if
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_area2, 1, int_var )
-  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2
+                     eh_area2(sn), 1, int_var )
+  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2(sn)
+
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     eh_area_tmp, 1, int_var )
+!  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area_tmp
+
   call CCTK_INFO(info_message)
 
+!  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_area2" )
+!  if ( ierr .lt. 0 ) then
+!    call CCTK_WARN ( 0, "cannot get lower bounds for area array" )
+!  end if
+!  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_area2" )
+!  if ( ierr .lt. 0 ) then
+!    call CCTK_WARN ( 0, "cannot get upper bounds for area array" )
+!  end if
+!  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_area2" )
+!  if ( ierr .lt. 0 ) then
+!    call CCTK_WARN ( 0, "cannot get local size for area array" )
+!  end if
+!
+!  if ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then
+!    eh_area2(sn-lbnd(1)) = eh_area_tmp
+!  end if
+!  print*,lbnd,ubnd,lsh
+
 end subroutine EHFinder_IntegrateArea
 
 
@@ -229,10 +254,15 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: int_var
+  CCTK_INT :: int_var, sn, k
   character(len=50) :: info_message
+!  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
+  CCTK_REAL :: centroid_x, centroid_y, centroid_z
 
   call CCTK_INFO ( 'Integrating centroid' )
+
+  sn = surface_counter - integrate_counter
+
   interp_x = center(1) + rsurf * sintheta * cosphi
   interp_y = center(2) + rsurf * sintheta * sinphi
   interp_z = center(3) + rsurf * costheta
@@ -244,31 +274,66 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
 
   int_tmp = interp_x * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid_x, 1, int_var )
+                     eh_centroid2_x(sn), 1, int_var )
 
   int_tmp = interp_y * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid_y, 1, int_var )
+                     eh_centroid2_y(sn), 1, int_var )
 
   int_tmp = interp_z * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid_z, 1, int_var )
-
-  eh_centroid_x = eh_centroid_x / eh_area2
-  eh_centroid_y = eh_centroid_y / eh_area2
-  eh_centroid_z = eh_centroid_z / eh_area2
-
-  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ',eh_centroid_x, &
-                                                             eh_centroid_y, &
-                                                             eh_centroid_z
+                     eh_centroid2_z(sn), 1, int_var )
+
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     centroid_x, 1, int_var )
+!
+!  int_tmp = interp_y * weights * da
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     centroid_y, 1, int_var )
+!
+!  int_tmp = interp_z * weights * da
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     centroid_z, 1, int_var )
+!
+!  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_area2" )
+!  if ( ierr .lt. 0 ) then
+!    call CCTK_WARN ( 0, "cannot get lower bounds for area array" )
+!  end if
+!  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_area2" )
+!  if ( ierr .lt. 0 ) then
+!    call CCTK_WARN ( 0, "cannot get upper bounds for area array" )
+!  end if
+!  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_area2" )
+!  if ( ierr .lt. 0 ) then
+!    call CCTK_WARN ( 0, "cannot get local size for area array" )
+!  end if
+
+    eh_centroid2_x(sn) = eh_centroid2_x(sn) / eh_area2(sn)
+    eh_centroid2_y(sn) = eh_centroid2_y(sn) / eh_area2(sn)
+    eh_centroid2_z(sn) = eh_centroid2_z(sn) / eh_area2(sn)
+
+!  if ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then
+!    k = sn - lbnd(1)
+!    eh_centroid2_x(k) = eh_centroid2_x(k) / eh_area2(k)
+!    eh_centroid2_y(k) = eh_centroid2_y(k) / eh_area2(k)
+!    eh_centroid2_z(k) = eh_centroid2_z(k) / eh_area2(k)
+!  end if
+
+  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', &
+                                       eh_centroid2_x(sn), &
+                                       eh_centroid2_y(sn), &
+                                       eh_centroid2_z(sn)
+
+!  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', &
+!                                       centroid_x, &
+!                                       centroid_y, &
+!                                       centroid_z
   call CCTK_INFO(info_message)
 
 end subroutine EHFinder_IntegrateCentroid
 
 
 ! This routine sets up the weights for the Simpsons rule integration
-
-! This routine sets up the weights for the Simpsons rule integration
 ! over the surface.
 subroutine EHFinder_InitWeights(CCTK_ARGUMENTS)
 
@@ -340,3 +405,39 @@ subroutine EHFinder_InitWeights(CCTK_ARGUMENTS)
   !  4/9  16/9  8/9  16/9  8/9  16/9  4/9
   !  1/9   4/9  2/9   4/9  2/9   4/9  1/9
 end subroutine EHFinder_InitWeights
+
+
+subroutine EHFinder_CopyArea(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  call CCTK_INFO('Copying area data')
+
+  eh_area = eh_area2
+
+end subroutine EHFinder_CopyArea
+
+
+subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  call CCTK_INFO('Copying centroid data')
+
+  eh_centroid_x = eh_centroid2_x
+  eh_centroid_y = eh_centroid2_y
+  eh_centroid_z = eh_centroid2_z
+
+end subroutine EHFinder_CopyCentroid