Second pass at finding points on the level set surface. It now uses

grid arrays and works in parallel. Still needs to be generalized to work
in other modes than full mode and an algorithm has to be implemented
to figure out how many surfaces are actually present in the level set
function.


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

4 files changed, 301 insertions, 19 deletions

diff --git a/interface.ccl b/interface.ccl
index ed19712..d1991ec 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -9,43 +9,59 @@ USES INCLUDE: MoL.h
 
 public:
 
-CCTK_REAL level_set type=GF TimeLevels=2
+CCTK_REAL level_set TYPE=GF TIMELEVELS=2
 {
   f
 } "The scalar level set function that defines the null surface"
 
-CCTK_REAL slevel_set type=GF TimeLevels=1
+CCTK_REAL slevel_set TYPE=GF TIMELEVELS=1
 {
   sf
 } "Source for the level set function"
 
-CCTK_REAL dlevel_set type=GF TimeLevels=1
+CCTK_REAL dlevel_set TYPE=GF TIMELEVELS=1
 {
   dfx, dfy, dfz, dfsq
 } "Derivatives of the level set function"
 
-CCTK_REAL ftmp_set type=GF TimeLevels=1
+CCTK_REAL ftmp_set TYPE=GF TIMELEVELS=1
 {
   ftmp, sftmp
 } "temporary variables used in pde re-parametrization"
 
-CCTK_REAL level_bak type=GF TimeLevels=1
+CCTK_REAL level_bak TYPE=GF TIMELEVELS=1
 {
   fbak
 } "Temporary placeholder for the level set during re-parametrization"
 
-CCTK_INT eh_mask_all type=GF TimeLevels=1
+CCTK_INT eh_mask_all TYPE=GF TIMELEVELS=1
 {
   eh_mask, tm_mask
 } "Masks to define active cells"
 
-CCTK_INT mask_bak type=GF TimeLevels=1
+CCTK_INT mask_bak TYPE=GF TIMELEVELS=1
 {
   eh_mask_bak
 } "Temporary placeholder for the mask during re-parametrization"
 
-CCTK_INT rep_mask type=GF TimeLevels=1
+CCTK_REAL surface_arrays TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=ntheta,nphi GHOSTSIZE=n_array_ghosts,n_array_ghosts DISTRIB=DEFAULT
+{
+  ctheta, cphi, rsurf, sinthetacosphi, sinthetasinphi, costheta
+} "Grid arrays for points on the surface"
+
+CCTK_REAL surface_tmp_arrays TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=ntheta,nphi GHOSTSIZE=n_array_ghosts,n_array_ghosts DISTRIB=DEFAULT
+{
+  drsurf, interp_x, interp_y, interp_z, f_interp,
+  dfdx_interp, dfdy_interp, dfdz_interp
+} "Tmeporary grid arrays for finding points on the surface"
+
+CCTK_INT surface_int_array TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=ntheta,nphi GHOSTSIZE=n_array_ghosts,n_array_ghosts DISTRIB=DEFAULT
+{
+  n_since_last_reduction
+} "Tmeporary integer grid array for finding points on the surface"
+
+CCTK_INT rep_mask TYPE=GF TIMELEVELS=1
 
-CCTK_INT re_param_control_pde type=SCALAR
-CCTK_INT re_param_control_approx type=SCALAR
-CCTK_REAL eh_area type=SCALAR
+CCTK_INT re_param_control_pde TYPE=SCALAR
+CCTK_INT re_param_control_approx TYPE=SCALAR
+CCTK_REAL eh_area TYPE=SCALAR
diff --git a/param.ccl b/param.ccl
index 16da9a2..60f1b2f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -30,22 +30,22 @@ KEYWORD initial_f "Initial surface choice"
 
 REAL initial_rad "Initial radius of surface"
 {
-(0.0:            :: "Positive please"
+0.0:            :: "Positive please"
 } 1.0
 
 REAL initial_a "Initial a coefficient of ellipsoid"
 {
-(0.0:            :: "Positive please"
+0.0:            :: "Positive please"
 } 1.0
 
 REAL initial_b "Initial b coefficient of ellipsoid"
 {
-(0.0:            :: "Positive please"
+0.0:            :: "Positive please"
 } 1.0
 
 REAL initial_c "Initial c coefficient of ellipsoid"
 {
-(0.0:            :: "Positive please"
+0.0:            :: "Positive please"
 } 1.0
 
 REAL rotation_alpha "Rotation angle around z-axis of ellipsoid"
@@ -90,7 +90,7 @@ REAL cas_b "Initial b coefficient of ovaloid of cassini"
 
 REAL shell_width "Width of the evolution region in units of the grid spacing"
 {
-(0.0:            :: "Positive please"
+0.0:            :: "Positive please"
 } 7.0
 
 BOOLEAN normalize "normalize the derivatives of the level set function"
@@ -133,7 +133,7 @@ KEYWORD re_param_int_method "Integration method in pde re-parametrization"
 
 INT re_param_max_iter "maximum number of iteration in the re-parametrization"
 {
-(0:        :: "Positive please"
+0:        :: "Positive please"
 } 800
 
 KEYWORD pde_differences "Type of finite diffencing used in pde re-parametrization"
@@ -155,7 +155,7 @@ INT reparametrize_every_approx "Re-parametrize every using approx method"
 
 INT last_iteration_number "The starting iteration number for the EH_Finder (last iteration number of the simulation)"
 {
-(0:              :: "Positive please"
+0:              :: "Positive please"
 } 0
 
 INT saved_iteration_every "How often was the numerical data saved"
@@ -163,6 +163,20 @@ INT saved_iteration_every "How often was the numerical data saved"
 1:               :: "Positive please"
 } 1
 
+INT ntheta "Number of points in the theta direction when finding points on the surface"
+{
+1:               :: "Positive and odd please"
+} 51
+
+INT nphi "Number of points in the phi direction when finding points on the surface"
+{
+1:               :: "Positive and odd please"
+} 51
+INT n_array_ghosts "Number of ghost points in the surface grid array"
+{
+1:               :: "Positive please"
+} 1
+
 shares: grid
 
 USES KEYWORD domain
diff --git a/schedule.ccl b/schedule.ccl
index 8cf7384..336f55b 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -15,6 +15,7 @@ else
   STORAGE: ftmp_set
   STORAGE: level_bak
   STORAGE: mask_bak
+  STORAGE: surface_arrays
 }
 
 # Check for metric_state
@@ -78,7 +79,7 @@ schedule EHFinder_Integrate at CCTK_ANALYSIS after EHFinder_FindSurface
 schedule EHFinder_FindSurface at CCTK_ANALYSIS
 {
   LANG: Fortran
-  STORAGE: eh_area
+  STORAGE: eh_area, surface_tmp_arrays, surface_int_array
   TRIGGER: eh_area
 } "Find Surface"
 
diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
new file mode 100644
index 0000000..e51542d
--- /dev/null
+++ b/src/EHFinder_FindSurface.F90
@@ -0,0 +1,251 @@
+! Finding the surface(s). Right now only on one processor and in full mode.
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k
+  CCTK_REAL, parameter :: eps = 1.0d-10
+  CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
+
+  CCTK_INT, dimension(4) :: bbox
+  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost
+
+  CCTK_REAL, dimension(3) :: center_loc, center
+  CCTK_INT :: nc_loc, nc
+  CCTK_REAL :: initial_radius, ncinv, maxdr_loc, maxdr, maxf_loc, maxf
+
+  CCTK_REAL :: dtheta, dphi, delta, rloc
+  CCTK_POINTER, dimension(3) :: interp_coords
+  CCTK_POINTER, dimension(4) :: out_array
+  CCTK_POINTER :: in_array
+  CCTK_INT, dimension(4), parameter :: op_indices = (/ 0, 0, 0, 0 /), &
+                                       op_codes = (/ 0, 1, 2, 3 /), &
+                                       out_types = (/ CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL /)
+  
+#include "include/physical_part.h"
+
+  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" )
+  end if
+!  print*,bbox
+  call CCTK_GroupgshGN ( status, cctkGH, 2, gsh, "ehfinder::surface_arrays" )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get global size for surface arrays" )
+  end if
+!  print*,gsh
+  call CCTK_GrouplbndGN ( status, cctkGH, 2, lbnd, "ehfinder::surface_arrays" )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get lower bounds for surface arrays" )
+  end if
+!  print*,lbnd
+  call CCTK_GroupubndGN ( status, cctkGH, 2, ubnd, "ehfinder::surface_arrays" )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get upper bounds for surface arrays" )
+  end if
+!  print*,ubnd
+  call CCTK_GrouplshGN ( status, cctkGH, 2, lsh, "ehfinder::surface_arrays" )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+  end if
+!  print*,lsh
+  call CCTK_GroupnghostzonesGN ( status, cctkGH, 2, nghost, "ehfinder::surface_arrays" )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+  end if
+  
+  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
+            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)
+            center_loc(3) = center_loc(3) + z(i,j,k)
+          end if
+        end if
+      end do
+    end do
+  end do
+
+  call CCTK_ReduceLocScalar ( status, cctkGH, -1, sum_handle, &
+                                    nc_loc, nc, CCTK_VARIABLE_INT )
+  call CCTK_ReduceLocArrayToArray1D ( status, cctkGH, -1, sum_handle, &
+                                      center_loc, center, 3, &
+                                      CCTK_VARIABLE_REAL )
+  ncinv = one / nc
+  center = ncinv * center
+  if ( symx ) center(1) = zero
+  if ( symy ) center(2) = zero
+  if ( symz ) center(3) = zero
+
+  dtheta = pi / ( ntheta - 1 )
+  dphi = two * pi / ( nphi - 1 )
+  do i = 1, lsh(1)
+    ctheta(i,:) = dtheta * ( i + lbnd(1) - 1 )
+  end do
+  do j = 1, lsh(2)
+    cphi(:,j) = dphi * ( j +lbnd(2) - 1 )
+  end do
+
+  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
+  
+  drsurf = two * delta
+  rsurf = zero
+  n_since_last_reduction = -1
+
+  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??" )
+  end if
+
+  call Util_TableCreateFromString ( table_handle, "order=3" )
+  if ( table_handle .lt. 0 ) then
+    call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
+  end if
+  
+  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
+
+  interp_coords(1) = CCTK_PointerTo(interp_x)
+  interp_coords(2) = CCTK_PointerTo(interp_y)
+  interp_coords(3) = CCTK_PointerTo(interp_z)
+
+  out_array(1) = CCTK_PointerTo(f_interp)
+
+  call CCTK_VarIndex ( in_array, "ehfinder::f" )
+
+! find starting point in all directions at the same time
+
+  find_starting_points: do
+
+    do j = 1, lsh(2)
+      do i = 1, lsh(1)
+        rloc = rsurf(i,j) + drsurf(i,j)
+        interp_x(i,j) = center(1) + rloc * sinthetacosphi(i,j)
+        interp_y(i,j) = center(2) + rloc * sinthetasinphi(i,j)
+        interp_z(i,j) = center(3) + rloc * costheta(i,j)
+        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
+      end do
+    end do
+
+    call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
+                                 interp_coords, 1, in_array, &
+                                 1, CCTK_VARIABLE_REAL, out_array(1) )
+                               
+    maxdr_loc = maxval ( drsurf )
+    call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
+                                    maxdr_loc, maxdr, CCTK_VARIABLE_REAL )
+    
+    if ( maxdr .lt. delta ) exit find_starting_points
+
+    do j = 1, lsh(2)
+      do i = 1, lsh(1)
+        if ( f_interp(i,j) .gt. zero ) then
+          drsurf(i,j) = half * drsurf(i,j)
+          n_since_last_reduction(i,j) = 0
+        else
+          rsurf(i,j) = rsurf(i,j) + drsurf(i,j)
+          if ( n_since_last_reduction(i,j) .ge. 1 ) then
+            drsurf(i,j) = half * drsurf(i,j)
+          end if
+        end if
+      end do
+    end do
+    
+  end do find_starting_points
+ 
+  out_array(2) = CCTK_PointerTo(dfdx_interp)
+  out_array(3) = CCTK_PointerTo(dfdy_interp)
+  out_array(4) = CCTK_PointerTo(dfdz_interp)
+
+  call Util_TableSetIntArray ( status, table_handle, 4, &
+                               op_indices, "operand_indices" )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+  endif
+  call Util_TableSetIntArray ( status, table_handle, 4, &
+                               op_codes, "operation_codes" )
+  if ( status .lt. 0 ) then
+    call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+  endif
+
+! find the surface points in all directions simultaneously.
+
+  find_surface_points: do
+
+    do j = 1, lsh(2)
+      do i = 1, lsh(1)
+        interp_x(i,j) = center(1) + rsurf(i,j) * sinthetacosphi(i,j)
+        interp_y(i,j) = center(2) + rsurf(i,j) * sinthetasinphi(i,j)
+        interp_z(i,j) = center(3) + rsurf(i,j) * costheta(i,j)
+      end do
+    end do
+
+    call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
+                                 interp_coords, 1, in_array, &
+                                 4, out_types, out_array )
+ 
+    maxf_loc = maxval ( abs ( f_interp ) )
+    call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
+                                    maxf_loc, maxf, CCTK_VARIABLE_REAL )
+    
+    if ( maxf .gt. eps ) then
+      do j = 1, lsh(2)
+        do i = 1, lsh(1)
+          rsurf(i,j) = rsurf(i,j) - f_interp(i,j) / &
+                       ( dfdx_interp(i,j) * sinthetacosphi(i,j) + &
+                         dfdy_interp(i,j) * sinthetasinphi(i,j) + &
+                         dfdz_interp(i,j) * costheta(i,j) )
+        end do
+      end do
+    else
+      exit find_surface_points
+    end if
+!    call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" )
+
+  end do find_surface_points
+
+!  print*,rsurf
+  call Util_TableDestroy ( status, table_handle )
+  
+end subroutine EHFinder_FindSurface