Put in some error checking so that the code doesn't keep working on calculating

surface quantities if finding points on the surface didn't work because of
interpolation problems.


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

4 files changed, 71 insertions, 17 deletions

diff --git a/interface.ccl b/interface.ccl
index 5890c5e..b061a46 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -49,6 +49,8 @@ CCTK_INT surface_index TYPE=GF TIMELEVELS=1
   sc
 }
 
+CCTK_INT find_surface_status TYPE=SCALAR
+
 CCTK_INT surface_integers TYPE=SCALAR
 {
   surface_counter
diff --git a/schedule.ccl b/schedule.ccl
index 38fec6b..3e5f227 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -19,6 +19,7 @@ else
   STORAGE: center_arrays
 #  STORAGE: eh_area, eh_centroid
   STORAGE: eh_area2, eh_centroid2, eh_circumference2
+  STORAGE: find_surface_status
 #  STORAGE: surface_index
 }
 
diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 8f1f1e7..bbd9521 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -46,8 +46,10 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 ! the 3D grid functions.
 #include "include/physical_part.h"
 
+! Initialize find_surface_status to 0. If something goes wrong it will be
+! set to -1.
+  find_surface_status = 0
   sn = surface_counter - integrate_counter + 1
-
   write(info_message,'(a26,i3)') 'Finding points on surface ',sn
   call CCTK_INFO ( info_message )
 
@@ -461,6 +463,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                  lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
                                  interp_coords, 1, in_array, &
                                  1, CCTK_VARIABLE_REAL, out_array(1) )
+    if ( status .lt. 0 ) then
+       call CCTK_INFO ( 'Interpolation failed. Giving up on finding surfaces' )
+       find_surface_status = -1
+       integrate_counter = integrate_counter - 1
+       return
+    end if
                                
     maxdr_loc = maxval ( drsurf )
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
@@ -526,6 +534,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                  lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
                                  interp_coords, 1, in_array, &
                                  4, out_types, out_array )
+    if ( status .lt. 0 ) then
+       call CCTK_INFO ( 'Interpolation failed. Giving up on finding surfaces' )
+       find_surface_status = -1
+       integrate_counter = integrate_counter - 1
+       return
+    end if
  
     maxf_loc = maxval ( abs ( f_interp ) )
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
index cc77255..81745ea 100644
--- a/src/EHFinder_Integrate2.F90
+++ b/src/EHFinder_Integrate2.F90
@@ -35,6 +35,11 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
                                                       CCTK_VARIABLE_REAL, &
                                                       CCTK_VARIABLE_REAL /)
   
+! If finding of surface failed don't try to integrate but exit.
+  if ( find_surface_status .lt. 0 ) then
+    return
+  endif
+
   call CCTK_INFO ( 'Finding surface element' )
   call CCTK_GroupbboxGN ( ierr, cctkGH, 4, bbox, "ehfinder::surface_arrays" )
   if ( ierr .lt. 0 ) then
@@ -65,12 +70,12 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
   if ( ierr .lt. 0 ) then
     call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
   end if
-  
+
   dtheta = ctheta(2,1) - ctheta(1,1)
   dphi = cphi(1,2) - cphi(1,1)
   dthetainv = one / dtheta
   dphiinv = one / dphi
-  
+
   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??" )
@@ -80,7 +85,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
   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 ??" )
@@ -207,6 +212,11 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
   character(len=30) :: info_message
   CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
 
+! If finding of surface failed don't try to integrate but exit.
+  if ( find_surface_status .lt. 0 ) then
+    return
+  endif
+
   call CCTK_INFO ( 'Integrating area' )
 
   sn = surface_counter - integrate_counter
@@ -264,6 +274,11 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
   CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
   CCTK_REAL :: centroid_x, centroid_y, centroid_z
 
+! If finding of surface failed don't try to integrate but exit.
+  if ( find_surface_status .lt. 0 ) then
+    return
+  endif
+
   call CCTK_INFO ( 'Integrating centroid' )
 
   sn = surface_counter - integrate_counter
@@ -367,6 +382,11 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS)
   CCTK_INT :: ithl, ithr, jphl, jphr, itha, jpha
   CCTK_REAL :: eq_circ_loc, eq_circ, pol_circ_loc, pol_circ
 
+! If finding of surface failed don't try to integrate but exit.
+  if ( find_surface_status .lt. 0 ) then
+    return
+  endif
+
   sn = surface_counter - integrate_counter
 
   ! Find out the lower bounds of the distributed integration grid arrays.
@@ -581,13 +601,18 @@ subroutine EHFinder_CopyArea(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  call CCTK_INFO('Copying area data')
-
 !  print*,Xeh_area0
 !  print*,Xeh_area20
 
-  eh_area = eh_area2
-
+! If finding of surface failed set area to zero.
+  print*,'Debug2 : ',find_surface_status
+  if ( find_surface_status .lt. 0 ) then
+    eh_area = zero
+    return
+  else
+    call CCTK_INFO('Copying area data')
+    eh_area = eh_area2
+  endif
 end subroutine EHFinder_CopyArea
 
 
@@ -601,14 +626,21 @@ subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  call CCTK_INFO('Copying centroid data')
-
 !  print*,Xeh_centroid0
 !  print*,Xeh_centroid20
 
-  eh_centroid_x = eh_centroid2_x
-  eh_centroid_y = eh_centroid2_y
-  eh_centroid_z = eh_centroid2_z
+! If finding of surface failed set centroids to zero.
+  if ( find_surface_status .lt. 0 ) then
+    eh_centroid_x = zero
+    eh_centroid_y = zero
+    eh_centroid_z = zero
+    return
+  else
+    call CCTK_INFO('Copying centroid data')
+    eh_centroid_x = eh_centroid2_x
+    eh_centroid_y = eh_centroid2_y
+    eh_centroid_z = eh_centroid2_z
+  end if
 
 end subroutine EHFinder_CopyCentroid
 
@@ -623,9 +655,14 @@ subroutine EHFinder_CopyCircumference(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  call CCTK_INFO('Copying circumference data')
-
-  eh_circ_eq = eh_circ_eq2
-  eh_circ_pol = eh_circ_pol2
+  if ( find_surface_status .lt. 0 ) then
+    eh_circ_eq = zero
+    eh_circ_pol = zero
+    return
+  else
+    call CCTK_INFO('Copying circumference data')
+    eh_circ_eq = eh_circ_eq2
+    eh_circ_pol = eh_circ_pol2
+  end if
 
 end subroutine EHFinder_CopyCircumference