Added more explanatory comments.

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

5 files changed, 142 insertions, 21 deletions

diff --git a/src/EHFinder_ReadData.F90 b/src/EHFinder_ReadData.F90
index 1076907..68516b6 100644
--- a/src/EHFinder_ReadData.F90
+++ b/src/EHFinder_ReadData.F90
@@ -1,10 +1,11 @@
-! Read metric, lapse and shift from files.
+! Read metric, lapse, shift and conformal factor from files.
 ! $Header$
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 
+!This routine reads in the metric (from ADMBase).
 subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -20,12 +21,25 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS)
   character(len=10) :: iteration_string
   CCTK_INT :: i, nc, ntot, res
 
+  ! Figure out which iteration number to read, based on the parameters
+  ! last_iteration_number (the last iteration in the numerical evolution
+  ! producing the metric) and saved_iteration_every (how often was the metric
+  ! saved) and the current iteration and save it in a string variable.
+
   write(iteration_string,'(i10)') last_iteration_number - &
                                   saved_iteration_every *cctk_iteration
+
+  ! Trim the string variable.
   iteration_string = adjustl(iteration_string)
   nc = len_trim(iteration_string)
 
+  ! Generate a string with the filenames. Note at present the requirement
+  ! therefore is that all iterations of one variable is written in the
+  ! same file.
   in_files = 'gxx gxy gxz gyy gyz gzz'
+
+  ! Fill in the string array, used to requesting the metric components at
+  ! the specified iteration number.
   var_names(1) = 'admbase::gxx[cctk_iteration='//iteration_string(1:nc)//']'
   var_names(2) = 'admbase::gxy[cctk_iteration='//iteration_string(1:nc)//']'
   var_names(3) = 'admbase::gxz[cctk_iteration='//iteration_string(1:nc)//']'
@@ -33,6 +47,7 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS)
   var_names(5) = 'admbase::gyz[cctk_iteration='//iteration_string(1:nc)//']'
   var_names(6) = 'admbase::gzz[cctk_iteration='//iteration_string(1:nc)//']'
 
+  ! merge all the varible names into a single string.
   in_vars = ' '
   ntot = 0
   do i = 1, 6
@@ -41,11 +56,14 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS)
     ntot = ntot + nc + 1
   end do
    
+  ! Call the routine that actualle does the file access. Note failures are
+  ! just silently ignored.
   call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars )
 
 end subroutine EHFinder_Read_Metric
 
 
+! This routine reads in the lapse (from ADMBase).
 subroutine EHFinder_Read_Lapse(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -73,6 +91,7 @@ subroutine EHFinder_Read_Lapse(CCTK_ARGUMENTS)
 end subroutine EHFinder_Read_Lapse
 
 
+! This routine reads in all the shift components (from ADMBAse).
 subroutine EHFinder_Read_Shift(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -110,6 +129,7 @@ subroutine EHFinder_Read_Shift(CCTK_ARGUMENTS)
 end subroutine EHFinder_Read_Shift
 
 
+! This routine reads in the conformal factor (from StaticConformal).
 subroutine EHFinder_Read_Conformal(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -137,6 +157,7 @@ subroutine EHFinder_Read_Conformal(CCTK_ARGUMENTS)
 end subroutine EHFinder_Read_Conformal
 
 
+! This routine reads in everything.
 subroutine EHFinder_ReadData(CCTK_ARGUMENTS)
 
   use EHFinder_mod
diff --git a/src/include/physical_part.h b/src/include/physical_part.h
index 6e396b9..40177a1 100644
--- a/src/include/physical_part.h
+++ b/src/include/physical_part.h
@@ -2,12 +2,14 @@
 ! $Header$
 
 ! nx, ny and nz contains the full size of the chunk.
+! nx, ny and nz are defined in EHFinder_mod.F90.
 
 nx = cctk_lsh(1)
 ny = cctk_lsh(2)
 nz = cctk_lsh(3)
 
 ! ngx, ngy, ngz contains the size of any ghostzones.
+! ngx, ngy and ngz are defined in EHFinder_mod.F90.
 
 ngx = cctk_nghostzones(1)
 ngy = cctk_nghostzones(2)
@@ -15,6 +17,7 @@ ngz = cctk_nghostzones(3)
 
 ! ixl and ixr is set to the minimum and maximum local cell index in the
 ! x-direction. The corresponding is done for the y- and z-directions.
+! ixl, ixr, jyl, jyr, kzl, kzr are defined in EHFinder_mod.F90.
 
 symx = .false.; symy = .false.; symz = .false.;
 
@@ -41,6 +44,8 @@ if ( cctk_bbox(6) .eq. 0 ) kzr = kzr - ngz
 ! boundaries. This is done by checking if the lowest index of the local grid
 ! as seen on the global grid is equal to zero. If this is the case adjust the
 ! minimum and maximum cell index appropriately with the size of the ghostzone.
+! Also symx, symy and symz are set to true.
+! symx, symy and symz are defined in EHFinder_mod.F90.
 
 if ( CCTK_EQUALS ( domain, 'octant' ) ) then
   if (cctk_lbnd(1) .eq. 0 ) then
@@ -115,3 +120,7 @@ if ( CCTK_EQUALS ( domain, 'bitant' ) ) then
     end if
   endif
 end if
+
+! At the end ixl, ixr, jyl, jyr, kzl, kzr should contain the minimum index 
+! and maximum index in each direction that are not ghost cells or
+! symmetry cells.
diff --git a/src/include/upwind_second.h b/src/include/upwind_second.h
index 8ab9ba8..f6004e9 100644
--- a/src/include/upwind_second.h
+++ b/src/include/upwind_second.h
@@ -1,44 +1,67 @@
-! Calculate second order upwinded differences.
+! Calculate second order upwinded differences and store them in the dfdx, dfdy
+! and dfdz grid functions.
 ! $Header$
 
+! First figure out the range of indices excluding ghost cells and symmetry
+! cells.
 # include "physical_part.h"
 
+! Set idx, idy and idz to 0.5*gridspacing
+! idx, idy and idz must be defined in the including routine.
 idx = half / cctk_delta_space(1)
 idy = half / cctk_delta_space(2)
 idz = half / cctk_delta_space(3)
 
-!if ( cctk_bbox(1) .eq. 0 ) eh_mask(1:ngx,:,:) = -2
-!if ( cctk_bbox(2) .eq. 0 ) eh_mask(nx-ngx+1:nx,:,:) = -2
-!if ( cctk_bbox(3) .eq. 0 ) eh_mask(:,1:ngy,:) = -2
-!if ( cctk_bbox(4) .eq. 0 ) eh_mask(:,ny-ngy+1:ny,:) = -2
-!if ( cctk_bbox(5) .eq. 0 ) eh_mask(:,:,1:ngz) = -2
-!if ( cctk_bbox(6) .eq. 0 ) eh_mask(:,:,nz-ngz+1:nz) = -2
-
+! Calculate appropriate one sided derivatives for all active cells.
+! Cells are active if eh_mask is greater or equal to zero. They are interior
+! cells if eh_mask == 0 and boundary cells if eh_mask > 0.
 do k = kzl, kzr
   do j = jyl, jyr
     do i = ixl, ixr
+
+      ! If this is an active cell...
       if ( eh_mask(i,j,k) .ge. 0 ) then
+
+        ! If this is not a boundary cell...
         if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
              ( .not. btest(eh_mask(i,j,k),1) ) ) then
+
+          ! If the cell to the left is not a boundary cell...
           if ( ( .not. btest(eh_mask(i-1,j,k),0) ) .or. &
                      ( eh_mask(i-1,j,k) .eq. -2 ) ) then
+
+            ! Calculate left handed one sided derivative
             al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
           else
             al = zero
           end if
+
+          ! If the cell to the right is not a boundary cell...
           if ( ( .not. btest(eh_mask(i+1,j,k),1) ) .or. &
                      ( eh_mask(i+1,j,k) .eq. -2 ) ) then 
+
+            ! Calculate right handed one sided derivative
             ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
           else
             ar = zero
           end if
+
+        ! Else if the cell is a left boundary cell.
         else if ( btest(eh_mask(i,j,k),0) ) then
+
+          ! calculate only right handed one sided difference.
           al = zero
           ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
+
+        ! Else it must be a right boundary cell.
         else
+
+          ! So calculate only the left handed one sided difference.
           al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
           ar = zero
         end if
+
+        ! Do the same for the y-derivative.
         if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
              ( .not. btest(eh_mask(i,j,k),3) ) ) then
           if ( ( .not. btest(eh_mask(i,j-1,k),2) ) .or. &
@@ -60,6 +83,8 @@ do k = kzl, kzr
           bl = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) )
           br = zero
         end if
+
+        ! And the z-derivative
         if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
              ( .not. btest(eh_mask(i,j,k),5) ) ) then
           if ( ( .not. btest(eh_mask(i,j,k-1),4) ) .or. &
@@ -82,21 +107,27 @@ do k = kzl, kzr
           cr = zero
         end if
 
+        ! Get the negative and positive parts of the one sided derivatives in
+        ! the x-direction.
         alminus = half*(al-abs(al))
         alplus = half*(al+abs(al))
         arminus = half*(ar-abs(ar))
         arplus = half*(ar+abs(ar))
 
+        ! Ditto for the y-direction.
         blminus = half*(bl-abs(bl))
         blplus = half*(bl+abs(bl))
         brminus = half*(br-abs(br))
         brplus = half*(br+abs(br))
 
+        ! Ditto for the z-direction.
         clminus = half*(cl-abs(cl))
         clplus = half*(cl+abs(cl))
         crminus = half*(cr-abs(cr))
         crplus = half*(cr+abs(cr))
 
+        ! If f>0 choose the correct one sided derivative depending on the
+        ! magnitude of the negative and positive parts
         if ( f(i,j,k) .gt. 0 ) then
           if ( abs(alplus) .gt. abs(arminus) ) then
             dfx(i,j,k) = alplus
@@ -115,6 +146,7 @@ do k = kzl, kzr
           end if
         endif
 
+        ! Ditto if f<0.
         if ( f(i,j,k) .lt. 0 ) then
           if ( abs(alminus) .gt. abs(arplus) ) then
             dfx(i,j,k) = alminus
@@ -132,23 +164,13 @@ do k = kzl, kzr
             dfz(i,j,k) = crplus
           end if
         end if
+
+      ! If the cell is not active set the derivatives to zero.
       else
         dfx(i,j,k) = zero
         dfy(i,j,k) = zero
         dfz(i,j,k) = zero
       end if
-!      if ( i .eq. 75 .and. j .eq. 40 .and. k .eq. 40 ) then
-!       print*,i,j,k
-!       print*,cl,cr
-!       print*,clminus,clplus
-!       print*,crminus,crplus
-!       print*,f(i,j,k)
-!       print*,dfz(i,j,k) 
-!       print*,eh_mask(i,j,k-2:k+2)
-!       print*,btest(eh_mask(i,j,k-2:k+2),4) 
-!       print*,btest(eh_mask(i,j,k-2:k+2),5) 
-!       pause
-!     end if
     end do
   end do
 end do
diff --git a/src/include/upwind_second2.h b/src/include/upwind_second2.h
index 6969eff..d2e3c17 100644
--- a/src/include/upwind_second2.h
+++ b/src/include/upwind_second2.h
@@ -1,28 +1,51 @@
 ! Calculate second order upwinded differences.
+! The only difference from upwind_second.h is that here the loop over cells
+! must be done in the including routine.
 ! $Header$
 
+! If this is an active cell...
 if ( eh_mask(i,j,k) .ge. 0 ) then
+
+  ! If this is not a boundary cell...
   if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
        ( .not. btest(eh_mask(i,j,k),1) ) ) then
+
+    ! If the cell to the left is not a boundary cell...
     if ( ( .not. btest(eh_mask(i-1,j,k),0) ) .or. &
                ( eh_mask(i-1,j,k) .eq. -2 ) ) then
+
+      ! Calculate left handed one sided derivative
       al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
     else
       al = zero
     end if
+
+    ! If the cell to the right is not a boundary cell...
     if ( ( .not. btest(eh_mask(i+1,j,k),1) ) .or. &
                ( eh_mask(i+1,j,k) .eq. -2 ) ) then 
+
+      ! Calculate right handed one sided derivative
       ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
     else
       ar = zero
     end if
+
+  ! Else if the cell is a left boundary cell.
   else if ( btest(eh_mask(i,j,k),0) ) then
+
+    ! calculate only right handed one sided difference.
     al = zero
     ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
+
+  ! Else it must be a right boundary cell.
   else
+
+    ! So calculate only the left handed one sided difference.
     al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
     ar = zero
   end if
+
+  ! Do the same for the y-derivative.
   if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
        ( .not. btest(eh_mask(i,j,k),3) ) ) then
     if ( ( .not. btest(eh_mask(i,j-1,k),2) ) .or. &
@@ -44,6 +67,8 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
     bl = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) )
     br = zero
   end if
+
+  ! And the z-derivative
   if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
        ( .not. btest(eh_mask(i,j,k),5) ) ) then
     if ( ( .not. btest(eh_mask(i,j,k-1),4) ) .or. &
@@ -66,21 +91,27 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
     cr = zero
   end if
 
+  ! Get the negative and positive parts of the one sided derivatives in
+  ! the x-direction.
   alminus = half*(al-abs(al))
   alplus = half*(al+abs(al))
   arminus = half*(ar-abs(ar))
   arplus = half*(ar+abs(ar))
 
+  ! Ditto for the y-direction.
   blminus = half*(bl-abs(bl))
   blplus = half*(bl+abs(bl))
   brminus = half*(br-abs(br))
   brplus = half*(br+abs(br))
 
+  ! Ditto for the z-direction.
   clminus = half*(cl-abs(cl))
   clplus = half*(cl+abs(cl))
   crminus = half*(cr-abs(cr))
   crplus = half*(cr+abs(cr))
 
+  ! If f>0 choose the correct one sided derivative depending on the
+  ! magnitude of the negative and positive parts
   if ( f(i,j,k) .gt. 0 ) then
     if ( abs(alplus) .gt. abs(arminus) ) then
       dfx(i,j,k) = alplus
@@ -99,6 +130,7 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
     end if
   endif
 
+  ! Ditto if f<0.
   if ( f(i,j,k) .lt. 0 ) then
     if ( abs(alminus) .gt. abs(arplus) ) then
       dfx(i,j,k) = alminus
@@ -116,6 +148,8 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
       dfz(i,j,k) = crplus
     end if
   end if
+
+! If the cell is not active set the derivatives to zero.
 else
   dfx(i,j,k) = zero
   dfy(i,j,k) = zero
diff --git a/src/include/upwind_shift_second2.h b/src/include/upwind_shift_second2.h
index 9965b0f..f5fff46 100644
--- a/src/include/upwind_shift_second2.h
+++ b/src/include/upwind_shift_second2.h
@@ -1,28 +1,59 @@
 ! Calculation of shift upwinded differences except at boundaries.
 ! $Header$
 
+! If this is an active cell...
 if ( eh_mask(i,j,k) .ge. 0 ) then
+
+  ! If the shift in the x-direction is <0.
   if ( betax(i,j,k) .lt. zero ) then
+
+    ! If both of the two nearest cells to the right are active cells...
     if ( ( eh_mask(i+1,j,k) .ge. 0 ) .and. ( eh_mask(i+2,j,k) .ge. 0 ) ) then
+
+      ! Calculate the right handed one sided derivative.
       dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
                            four * f(i+1,j,k) - f(i+2,j,k) )
+
+    ! Else if only the nearest neighbour to the right is active...
     else if ( eh_mask(i+1,j,k) .ge. 0 ) then
+
+      ! Calculate the centered derivative.
       dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
+
+    ! Else it must be a boundary cell...
     else
+
+      ! So calculate the left handed one sided derivative.
       dfx(i,j,k) = idx * ( three * f(i,j,k) - &
                            four * f(i-1,j,k) + f(i-2,j,k) )
     end if
+
+  ! Else if the shift is >= 0.
   else
+
+    ! If both of the two nearest cells to the left are active cells...
     if ( ( eh_mask(i-1,j,k) .ge. 0 ) .and. ( eh_mask(i-2,j,k) .ge. 0 ) ) then
+
+      ! Calculate the left handed one sided derivative.
       dfx(i,j,k) = idx * ( three * f(i,j,k) - &
                            four * f(i-1,j,k) + f(i-2,j,k) )
+
+    ! Else if only the nearest neighbour to the left is active...
     else if ( eh_mask(i-1,j,k) .ge. 0 ) then
+
+      ! Calculate the centered derivative.
       dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
+
+    ! Else it must be a boundary cell...
     else
+
+      ! So calculate the left handed one sided derivative.
       dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
                            four * f(i+1,j,k) - f(i+2,j,k) )
     end if
   end if
+
+  ! Ditto for the y-derivative.
   if ( betay(i,j,k) .lt. zero ) then
     if ( ( eh_mask(i,j+1,k) .ge. 0 ) .and. ( eh_mask(i,j+2,k) .ge. 0 ) ) then
       dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
@@ -44,6 +75,8 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
                            four * f(i,j+1,k) - f(i,j+2,k) )
     end if
   end if
+
+  ! Ditto for the z-derivative.
   if ( betaz(i,j,k) .lt. zero ) then
     if ( ( eh_mask(i,j,k+1) .ge. 0 ) .and. ( eh_mask(i,j,k+2) .ge. 0 ) ) then
       dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
@@ -65,6 +98,8 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
                            four * f(i,j,k+1) - f(i,j,k+2) )
     end if
   end if
+
+! If the cell is not active set the derivatives to zero.
 else
   dfx(i,j,k) = zero
   dfy(i,j,k) = zero