Major modification of almost all source files to use vector groups for various

variables. Not tested in all detail, but the standard features seem to work.
The version before all these changes was tagged with PRE_MULTI.


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

1 files changed, 62 insertions, 99 deletions

diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 95e0e63..2537408 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -15,15 +15,16 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: i, j, k, im, jm, sn
+  CCTK_INT :: i, j, k, l, 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
+  character(len=80) :: info_message
   character(len=200) :: surface_interp
   CCTK_INT :: surface_interp_len
   character(len=7) :: surface_order
+  character(len=15) :: vname
 
   CCTK_INT, dimension(4) :: bbox
   CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost, maxl
@@ -54,7 +55,9 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 ! set to -1.
   find_surface_status = 0
   sn = surface_counter - integrate_counter + 1
-  write(info_message,'(a26,i3)') 'Finding points on surface ',sn
+  l = levelset_counter
+  write(info_message,'(a26,i3,a14,i3)') 'Finding points on surface ', sn, &
+                                 ' in level set ',l
   call CCTK_INFO ( info_message )
 
   ! Find information about the local and global properties of the 2D
@@ -103,13 +106,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     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 ) ) .and. &
+          if ( eh_mask(i,j,k,l) .eq. 0 ) then
+            if ( ( ( f(i,j,k,l) * f(i-1,j,k,l) .lt. zero ) .or. &
+                   ( f(i,j,k,l) * f(i+1,j,k,l) .lt. zero ) .or. &
+                   ( f(i,j,k,l) * f(i,j-1,k,l) .lt. zero ) .or. &
+                   ( f(i,j,k,l) * f(i,j+1,k,l) .lt. zero ) .or. &
+                   ( f(i,j,k,l) * f(i,j,k-1,l) .lt. zero ) .or. &
+                   ( f(i,j,k,l) * f(i,j,k+1,l) .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. &
@@ -352,16 +355,6 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   theta_sym_factor = two * pi / ( utheta - ltheta )
   phi_sym_factor = two * pi / ( uphi - lphi )
 
-!  print*
-!  print*,ltheta/pi,utheta/pi,lphi/pi,uphi/pi
-!  print*,sym_factor
-!  print*,center
-!  print*,symx,symy,symz
-!  ! For now symmetries are not handled properly
-!  if ( symx ) center(1) = zero
-!  if ( symy ) center(2) = zero
-!  if ( symz ) center(3) = zero
-
   ! 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.
@@ -384,13 +377,6 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   costheta = cos(ctheta)
   sinphi = sin(cphi)
   cosphi = cos(cphi)
-!  do j = 1, lsh(2)
-!    do i = 1, lsh(1)
-!      sinthetacosphi(i,j) = sin(ctheta(i,j)) * cos(cphi(i,j))
-!      sinthetasinphi(i,j) = sin(ctheta(i,j)) * sin(cphi(i,j))
-!      costheta(i,j) = cos(ctheta(i,j))
-!    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
@@ -442,7 +428,8 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   out_array(1) = CCTK_PointerTo(f_interp)
 
   ! Get the variable index for f
-  call CCTK_VarIndex ( in_array, "ehfinder::f" )
+  write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1, ']'
+  call CCTK_VarIndex ( in_array, vname )
 
 ! Loop to find starting point in all directions at the same time. For
 ! simplicity all ghost zones are found on all processors
@@ -492,9 +479,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                     maxf_loc, maxf, CCTK_VARIABLE_REAL )
     
     maxl = maxloc ( abs ( f_interp ) )
-!    print*,maxl,rsurf(maxl(1),maxl(2)),&
-!           drsurf(maxl(1),maxl(2)),&
-!           f_interp(maxl(1),maxl(2))
+
     if ( maxdr .lt. min_delta ) exit find_starting_points
 
     do j = 1, lsh(2)
@@ -559,14 +544,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                     maxf_loc, maxf, CCTK_VARIABLE_REAL )
     
     maxl = maxloc ( abs ( f_interp ) )
-!    print*,k,maxf
-!    print*,maxf,maxl,f_interp(maxl(1),maxl(2)),rsurf(maxl(1),maxl(2)), &
-!                     dfdx_interp(maxl(1),maxl(2))*sintheta(maxl(1),maxl(2))* &
-!                                                  cosphi(maxl(1),maxl(2))+ &
-!                     dfdy_interp(maxl(1),maxl(2))*sintheta(maxl(1),maxl(2))* &
-!                                                  sinphi(maxl(1),maxl(2))+ &
-!                     dfdz_interp(maxl(1),maxl(2))*costheta(maxl(1),maxl(2))
-!    print*
+
     if ( maxf .gt. eps ) then
       do j = 1, lsh(2)
         do i = 1, lsh(1)
@@ -694,6 +672,42 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 end subroutine EHFinder_FindSurface
 
 
+subroutine EHFinder_LevelSetLoopInit(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  more_levelsets = 1
+  levelset_counter = 1
+
+end subroutine EHFinder_LevelSetLoopInit
+
+
+subroutine EHFinder_UpdateCounter(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  if ( levelset_counter .lt. eh_number_level_sets ) then
+    levelset_counter = levelset_counter + 1
+    more_levelsets = 1
+  else
+    more_levelsets = 0
+  end if
+
+end subroutine EHFinder_UpdateCounter
+
+
 subroutine EHFinder_CountSurfacesInit(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -732,12 +746,12 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
 
   ! Find the location of the minimum among active points not already found
   ! to be inside a surface.
-  minl = minloc ( f(ixl:ixr,jyl:jyr,kzl:kzr), &
+  minl = minloc ( f(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter), &
                ( sc(ixl:ixr,jyl:jyr,kzl:kzr) .eq. 0 ) .and. &
-               ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 ) )
+               ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) )
 
   ! Find the value of f at that location.
-  minf_loc = f(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1)
+  minf_loc = f(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1,levelset_counter)
 
   ! Find the minimum over all processors.
   call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, &
@@ -810,8 +824,8 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
     do j = jyl, jyr
       do i = ixl, ixr
         marked = .false.
-        if ( ( eh_mask(i,j,k) .ge. 0 ) .and. &
-             ( f(i,j,k) .lt. 0 ) .and. &
+        if ( ( eh_mask(i,j,k,levelset_counter) .ge. 0 ) .and. &
+             ( f(i,j,k,levelset_counter) .lt. 0 ) .and. &
              ( sc(i,j,k) .ne. surface_counter ) ) then
 
           ! Find out if we are on the boundary
@@ -857,58 +871,6 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
             marked = .true.
           end if
 
-!          if ( any ( sc(i+i1:i+i2,j+j1:j+j2,k+k1:k+k2) &
-!                                       .eq. surface_counter ) ) then
-!            marked = .true.
-!          end if
-
-!          ! First look in the -x direction.
-!          if ( .not. btest(eh_mask(i,j,k),0) ) then
-!            if ( sc(i-1,j,k) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!            if ( .not. btest(eh_mask(i-1,j,k),2) ) then
-!              if ( sc(i-1,j-1,k) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!            if ( .not. btest(eh_mask(i-1,j,k),3) ) then
-!              if ( sc(i-1,j+1,k) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!            if ( .not. btest(eh_mask(i-1,j,k),4) ) then
-!              if ( sc(i-1,j,k-1) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!            if ( .not. btest(eh_mask(i-1,j,k),5) ) then
-!              if ( sc(i-1,j,k+1) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!          end if
-!          if ( .not. btest(eh_mask(i,j,k),1) ) then
-!            if ( sc(i+1,j,k) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!          end if
-!          if ( .not. btest(eh_mask(i,j,k),2) ) then
-!            if ( sc(i,j-1,k) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!          end if
-!          if ( .not. btest(eh_mask(i,j,k),3) ) then
-!            if ( sc(i,j+1,k) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!          end if
-!          if ( .not. btest(eh_mask(i,j,k),4) ) then
-!            if ( sc(i,j,k-1) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!          end if
-!          if ( .not. btest(eh_mask(i,j,k),5) ) then
-!            if ( sc(i,j,k+1) .eq. surface_counter ) then
-!              marked = .true.
-!            end if
-!          end if
           if ( marked ) then
             sc(i,j,k) = surface_counter
             n_loc = n_loc + 1
@@ -958,13 +920,14 @@ subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS)
 
   ! Write out how many surfaces were found.
   if ( surface_counter .eq. 1 ) then
-    write(info_message,'(a10,i3,a8)') 'There are ',surface_counter,' surface'
+    write(info_message,'(a9,i3,a22,i3)') 'There is ',surface_counter, &
+                                      ' surface in level set ',levelset_counter
   else
-    write(info_message,'(a10,i3,a9)') 'There are ',surface_counter,' surfaces'
+    write(info_message,'(a10,i3,a23,i3)') 'There are ',surface_counter, &
+                                      ' surfaces in level set ',levelset_counter
   end if
   call CCTK_INFO(info_message)
 
   integrate_counter = surface_counter
 
-!  print*,'EHFinder_InfoSurfaces : ',integrate_counter
 end subroutine EHFinder_InfoSurfaces