diff options
| -rw-r--r-- | interface.ccl | 2 | ||||
| -rw-r--r-- | schedule.ccl | 1 | ||||
| -rw-r--r-- | src/EHFinder_FindSurface.F90 | 16 | ||||
| -rw-r--r-- | src/EHFinder_Integrate2.F90 | 69 |
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 |