Added more comments.

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

1 files changed, 45 insertions, 3 deletions

diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 4a62ae0..6a2e401 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -1,4 +1,4 @@
-! Finding the surface(s). Right now only in full mode.
+! Finding the surface(s). Right now only in full mode but fully parallelised.
 ! $Header$
 
 #include "cctk.h"
@@ -38,8 +38,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                                       CCTK_VARIABLE_REAL, &
                                                       CCTK_VARIABLE_REAL /)
   
+! Find index ranges for points excluding ghost zones and symmetry points for
+! the 3D grid functions.
 #include "include/physical_part.h"
 
+
+  ! Find information about the local and global properties of the 2D
+  ! Grid arrays.
   call CCTK_GroupbboxGN ( status, cctkGH, 4, bbox, "ehfinder::surface_arrays" )
   if ( status .lt. 0 ) then
     call CCTK_WARN ( 0, "cannot get bounding box for surface arrays" )
@@ -70,11 +75,15 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
   end if
   
+  ! If the user defined center point should be used...
   if ( use_user_center ) 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.
     nc_loc = 0; center_loc = zero
     delta = minval ( cctk_delta_space )
     do k = kzl, kzr
@@ -97,6 +106,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       end do
     end do
 
+    ! Reduce over the processors.
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, sum_handle, &
                                     nc_loc, nc, CCTK_VARIABLE_INT )
     call CCTK_ReduceLocArrayToArray1D ( status, cctkGH, -1, sum_handle, &
@@ -105,6 +115,8 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     ncinv = one / nc
     center = ncinv * center
 !    print*,center
+
+    ! For now symmetries are not handled properly
     if ( symx ) center(1) = zero
     if ( symy ) center(2) = zero
     if ( symz ) center(3) = zero
@@ -112,6 +124,9 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 !  print*,'Center : ',center
 !  print*,'use_user_center : ',use_user_center
 
+  ! Find dtheta and dphi and initialise the theta and phi grid arrays.
+  ! Here i + lbnd(1) - 1 is the global index for theta and
+  !      j + lbnd(2) - 1 is the global index for phi.
   dtheta = pi / ( ntheta - 1 )
   dphi = two * pi / ( nphi - 1 )
   dthetainv = one / dtheta
@@ -123,6 +138,8 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     cphi(:,j) = dphi * ( j +lbnd(2) - 1 )
   end do
 
+  ! Calculate sines and cosines and store them in grid arrays since they
+  ! are expensive.
   sintheta = sin(ctheta)
   costheta = cos(ctheta)
   sinphi = sin(cphi)
@@ -135,10 +152,19 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 !    end do
 !  end do
   
+  ! rsurf is the radial grid array (as function of theta and phi) and is
+  ! initialised to zero. The grid array drsurf contains the changes to rsurf
+  ! and is initialised to 2*Delta. The grid array n_since_last_reduction is
+  ! a help in keeping track on how long ago drsurf for a given point has been
+  ! reduced and is initialised to -1 to indicate that no reduction has taken
+  ! place.
   drsurf = two * delta
   rsurf = zero
   n_since_last_reduction = -1
 
+  ! Get and interpolation handle for Hermite polynomial interpolation. This
+  ! is used to avoid problems with non-continuous derivatives in the Newton
+  ! iteration later on.
 !  call CCTK_InterpHandle ( interp_handle, "Lagrange polynomial interpolation" )
   call CCTK_InterpHandle ( interp_handle, "Hermite polynomial interpolation" )
   if ( interp_handle .lt. 0 ) then
@@ -146,36 +172,52 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
  &  Forgot to activate an implementation providing interpolation operators??" )
   end if
 
+  ! For now order=2 is hard wired.
 !  call Util_TableCreateFromString ( table_handle, "order=3" )
   call Util_TableCreateFromString ( table_handle, "order=2" )
   if ( table_handle .lt. 0 ) then
     call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
   end if
   
+  ! 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 ??" )
   endif
 
+  ! Get pointers to the grid arrays containing the interpolation points.
   interp_coords(1) = CCTK_PointerTo(interp_x)
   interp_coords(2) = CCTK_PointerTo(interp_y)
   interp_coords(3) = CCTK_PointerTo(interp_z)
 
+  ! Get pointer to the interpolation output array
   out_array(1) = CCTK_PointerTo(f_interp)
 
+  ! Get the variable index for f
   call CCTK_VarIndex ( in_array, "ehfinder::f" )
 
-! find starting point in all directions at the same time
-
+! Loop to find starting point in all directions at the same time. For
+! simplicity all ghost zones are found on all processors
+! The algorithm simply starts at the center and works it way outwards
+! in all directions. When the function changes signs, it reduces the
+! changes in radius until all points have negative f values and are
+! within one grid cells difference from the surface.
   find_starting_points: do
 
     do j = 1, lsh(2)
       do i = 1, lsh(1)
+
+        ! The radius that are tested in this direction.
         rloc = rsurf(i,j) + drsurf(i,j)
+
+        ! Calculate the interpolation points.
         interp_x(i,j) = center(1) + rloc * sintheta(i,j) * cosphi(i,j)
         interp_y(i,j) = center(2) + rloc * sintheta(i,j) * sinphi(i,j)
         interp_z(i,j) = center(3) + rloc * costheta(i,j)
+
+        ! If the point has been reduced once, increment n_since_last_reduction.
+        ! 
         if ( n_since_last_reduction(i,j) .ge. 0 ) then 
           n_since_last_reduction(i,j) = n_since_last_reduction(i,j) + 1
         end if