Added more comments.

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

2 files changed, 83 insertions, 9 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
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index 4aa43ff..c8e2f19 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -1,4 +1,4 @@
-! Initialisation of the level set function
+! Initialisation of the level set function and various other things.
 ! $Header$
 
 #include "cctk.h"
@@ -21,31 +21,51 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
   CCTK_REAL, dimension(3,3) :: txyz
   CCTK_REAL :: cosa, sina, cosb, sinb, cosc, sinc
 
+  ! Set up the value used in interiour inactive cells.
   ex_value = - ( one + shell_width ) * maxval(cctk_delta_space)
 
+  ! Initialize the current_iteration variable.
   current_iteration = last_iteration_number
 
+  ! Initialize the last_time varible. Note the parameters should be
+  ! chosen to be consistent with the run producing the numerical data.
+  ! F.ex. if dt was 0.1 but the data was only stored every 4 iterations, the
+  ! parameters should be chosen so that dt now is 0.4.
   last_time = abs(cctk_delta_time) * last_iteration_number / &
                                      saved_iteration_every
-  
+
   cctk_time = last_time
 
+  ! Initialize nx, ny and nz based on the local gris size. They are defined
+  ! in the module EHFinder_mod and will therefore be known from now on by all
+  ! routine that uses this module.
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
   nz = cctk_lsh(3)
 
+  ! If a sphere is requested...
   if ( CCTK_EQUALS( initial_f, 'sphere' ) ) then
+
+    ! Set up a sphere of radius initial_rad and translated 
+    ! (translate_x,translate_y,translate_z) away from the origin.
     f = sqrt( ( x - translate_x )**2 + &
               ( y - translate_y )**2 + &
               ( z - translate_z )**2 ) - initial_rad
   end if
+
+  ! If an ellipsoid is requested...
   if ( CCTK_EQUALS( initial_f, 'ellipsoid' ) ) then
+
+    ! Calculate sines and cosines of the rotation parameters.
     cosa = cos(rotation_alpha)
     sina = sin(rotation_alpha)
     cosb = cos(rotation_beta)
     sinb = sin(rotation_beta)
     cosc = cos(rotation_gamma)
     sinc = sin(rotation_gamma)
+
+    ! Set up the rotation matrix. The order is alpha around the z-axis,
+    ! beta around the y-axis and finally gamma around the x-axis.
     txyz(1,1) = cosa * cosb
     txyz(1,2) = sina * cosb
     txyz(1,3) = -sinb
@@ -55,7 +75,11 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
     txyz(3,1) = cosa * sinb * cosc + sina * sinc
     txyz(3,2) = sina * sinb * cosc - cosa * sinc
     txyz(3,3) = cosb * cosc
-    print*,txyz
+!    print*,txyz
+
+    ! Apply the rotations and translation for all points on the grid.
+    ! Even though at first glance it looks like the translation is done
+    ! first, the opposite is actually true.
     do k = 1, nz
       do j = 1, ny
         do i = 1, nx
@@ -70,18 +94,26 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
       end do
     end do
   end if
+
+  ! if an ovaloid of Cassini is requested...
   if ( CCTK_EQUALS( initial_f, 'cassini' ) ) then
     f = (x**2+y**2+z**2)**2 + cas_a**4 - &
         2*cas_a**2*(x**2 - (y**2+z**2)) - cas_b**4
   end if
 
+  ! Initialise the internal mask.
   eh_mask = 0
+
+  ! Find the maximal grid spacing.
   delta = maxval ( cctk_delta_space )
 
+  ! Initialise the ghostzone information. These variables are defined in
+  ! the module EHFinder_mod.
   ngx = cctk_nghostzones(1)
   ngy = cctk_nghostzones(2)
   ngz = cctk_nghostzones(3)
 
+  ! Get handles for various reduction operations.
   call CCTK_ReductionArrayHandle ( max_handle, 'maximum' )
   if ( max_handle .lt. 0 ) then
     call CCTK_WARN(0,'Could not obtain a handle for maximum reduction')
@@ -95,7 +127,9 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
     call CCTK_WARN(0,'Could not obtain a handle for sum reduction')
   end if
 
-
+  ! Register a coodinatsystem for the 2D surface grid arrays. This has no
+  ! effect yet, but should do something useful, when TAGS can be used to
+  ! assign coordinates to grid arrays.
   call CCTK_CoordRegisterSystem ( ierr, 2, 'cart2d' )
   if ( ierr .lt. 0 ) then
     call CCTK_WARN(0,'Could not register a 2D coordinate system as "cart2d"')
@@ -113,7 +147,5 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
 
   call CCTK_INFO('2d coordinate system registered')
 
-! set up weights for Simpsons rule integration
-
   return
 end subroutine EHFinder_Init