! Finding the surface(s). Right now only in full mode but fully parallelised. ! $Header$ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) use EHFinder_mod use EHFinder_fuzzy implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS 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=80) :: info_message character(len=128) :: warn_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 CCTK_REAL, dimension(3) :: center_loc CCTK_INT :: nc_loc, nc CCTK_REAL :: initial_radius, ncinv, maxdr_loc, maxdr, maxf_loc, maxf CCTK_REAL, dimension(3) :: crange_min_loc, crange_max_loc, & crange_min_glob, crange_max_glob CCTK_REAL :: ltheta, utheta, lphi, uphi CCTK_REAL :: dtheta, dphi, dthetainv, dphiinv, min_delta, rloc CCTK_POINTER, dimension(3) :: interp_coords CCTK_POINTER, dimension(4) :: out_array CCTK_INT :: 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 /) ! Find index ranges for points excluding ghost zones and symmetry points for ! 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 ! Set the index to the current level set function. l = levelset_counter ! Write an info message about what is going on. write(info_message,'(a26,i3,a14,i3)') 'Finding points on surface ', sn, & ' in level set ',l call CCTK_INFO ( trim(info_message) ) ! Find information about the local and global properties of the 2D ! Grid arrays. 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 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 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 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 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 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 ! If the user defined center point should be used... if ( use_user_center .gt. 0 ) then center(1) = center_x center(2) = center_y center(3) = center_z else ! If not find an approximation for it by averaging the coodinates that ! are within one grid cell distance from the surface in question. nc_loc = 0; center_loc = zero min_delta = minval ( cctk_delta_space ) do k = kzl, kzr do j = jyl, jyr do i = ixl, ixr 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. & ( fuzzy_eq ( sc(i,j,k), real(sn,wp), small ) .or. & fuzzy_eq ( sc(i-1,j,k), real(sn,wp), small ) .or. & fuzzy_eq ( sc(i+1,j,k), real(sn,wp), small ) .or. & fuzzy_eq ( sc(i,j-1,k), real(sn,wp), small ) .or. & fuzzy_eq ( sc(i,j+1,k), real(sn,wp), small ) .or. & fuzzy_eq ( sc(i,j,k-1), real(sn,wp), small ) .or. & fuzzy_eq ( sc(i,j,k+1), real(sn,wp), small ) ) ) 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 ! Reduce over the processors. 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 end if ! Try to figure out the symmetries. First find the range in coordinates ! for the points inside the current surface. crange_min_loc(1) = minval ( x, mask = fuzzy_eq ( sc, real(sn,wp), small ) ) crange_min_loc(2) = minval ( y, mask = fuzzy_eq ( sc, real(sn,wp), small ) ) crange_min_loc(3) = minval ( z, mask = fuzzy_eq ( sc, real(sn,wp), small ) ) crange_max_loc(1) = maxval ( x, mask = fuzzy_eq ( sc, real(sn,wp), small ) ) crange_max_loc(2) = maxval ( y, mask = fuzzy_eq ( sc, real(sn,wp), small ) ) crange_max_loc(3) = maxval ( z, mask = fuzzy_eq ( sc, real(sn,wp), small ) ) ! Reduce over all processors. call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & crange_min_loc, crange_min_glob, & 3, CCTK_VARIABLE_REAL ) call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, & crange_max_loc, crange_max_glob, & 3, CCTK_VARIABLE_REAL ) ! Full mode. Modify if necessary. ltheta = zero; utheta = pi lphi = zero; uphi = two * pi sym_factor = one s_symx = .false.; s_symy = .false.; s_symz = .false. if ( CCTK_EQUALS ( domain, 'bitant' ) ) then ! Bitant with symmetry in the z-direction if ( CCTK_EQUALS ( bitant_plane, 'xy' ) ) then ! If the minimum z-coordinate of the points inside the surface is less ! than zero, then the surface must share the grid symmetry across the ! xy-plane. if ( crange_min_glob(3) .lt. zero ) then ltheta = zero; utheta = half * pi lphi = zero; uphi = two * pi sym_factor = two center(3) = zero s_symz = .true. end if end if ! Bitant with symmetry in the y-direction if ( CCTK_EQUALS ( bitant_plane, 'xz' ) ) then if ( crange_min_glob(2) .lt. zero ) then ltheta = zero; utheta = pi lphi = zero; uphi = pi sym_factor = two center(2) = zero s_symy = .true. end if end if ! Bitant with symmetry in the x-direction if ( CCTK_EQUALS ( bitant_plane, 'yz' ) ) then if ( crange_min_glob(1) .lt. zero ) then ltheta = zero; utheta = pi lphi = -half * pi; uphi = half * pi sym_factor = two center(1) = zero s_symx = .true. end if end if end if if ( CCTK_EQUALS ( domain, 'quadrant' ) ) then ! Quadrant in x-direction if ( CCTK_EQUALS ( quadrant_direction, 'x' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & ( crange_min_glob(2) .lt. zero ) ) then ltheta = zero; utheta = half * pi lphi = zero; uphi = pi sym_factor = four center(3) = zero; center(2) = zero s_symz = .true.; s_symy = .true. else if ( crange_min_glob(3) .lt. zero ) then ltheta = zero; utheta = half * pi lphi = zero; uphi = two * pi sym_factor = two center(3) = zero s_symz = .true. else if ( crange_min_glob(2) .lt. zero ) then ltheta = zero; utheta = pi lphi = zero; uphi = pi sym_factor = two center(2) = zero s_symy = .true. end if end if ! Quadrant in y-direction if ( CCTK_EQUALS ( quadrant_direction, 'y' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & ( crange_min_glob(1) .lt. zero ) ) then ltheta = zero; utheta = half * pi lphi = -half * pi; uphi = half * pi sym_factor = four center(3) = zero; center(1) = zero s_symz = .true.; s_symx = .true. else if ( crange_min_glob(3) .lt. zero ) then ltheta = zero; utheta = half * pi lphi = zero; uphi = two * pi sym_factor = two center(3) = zero s_symz = .true. else if ( crange_min_glob(1) .lt. zero ) then ltheta = zero; utheta = pi lphi = -half * pi; uphi = half * pi sym_factor = two center(1) = zero s_symx = .true. end if end if ! Quadrant in z-direction if ( CCTK_EQUALS ( quadrant_direction,'z' ) ) then if ( ( crange_min_glob(2) .lt. zero ) .and. & ( crange_min_glob(1) .lt. zero ) ) then ltheta = zero; utheta = pi lphi = zero; uphi = half * pi sym_factor = four center(2) = zero; center(1) = zero s_symy = .true.; s_symx = .true. else if ( crange_min_glob(2) .lt. zero ) then ltheta = zero; utheta = pi lphi = zero; uphi = pi sym_factor = two center(2) = zero s_symy = .true. else if ( crange_min_glob(1) .lt. zero ) then ltheta = zero; utheta = pi lphi = -half * pi; uphi = half * pi sym_factor = two center(1) = zero s_symx = .true. end if end if end if ! Octant if ( CCTK_EQUALS ( domain, 'octant' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & ( crange_min_glob(2) .lt. zero ) .and. & ( crange_min_glob(1) .lt. zero ) ) then ltheta = zero; utheta = half * pi lphi = zero; uphi = half * pi sym_factor = eight center(3) = zero; center(2) = zero; center(1) = zero s_symz = .true.; s_symy = .true.; s_symx = .true. else if ( ( crange_min_glob(3) .lt. zero ) .and. & ( crange_min_glob(2) .lt. zero ) ) then ltheta = zero; utheta = half * pi lphi = zero; uphi = pi sym_factor = four center(3) = zero; center(2) = zero s_symz = .true.; s_symy = .true. else if ( ( crange_min_glob(3) .lt. zero ) .and. & ( crange_min_glob(1) .lt. zero ) ) then ltheta = zero; utheta = half * pi lphi = -half * pi; uphi = half * pi sym_factor = four center(3) = zero; center(1) = zero s_symz = .true.; s_symx = .true. else if ( ( crange_min_glob(2) .lt. zero ) .and. & ( crange_min_glob(1) .lt. zero ) ) then ltheta = zero; utheta = pi lphi = zero; uphi = half * pi sym_factor = four center(2) = zero; center(1) = zero s_symy = .true.; s_symx = .true. else if ( crange_min_glob(3) .lt. zero ) then ltheta = zero; utheta = half * pi lphi = zero; uphi = two * pi sym_factor = two center(3) = zero s_symz = .true. else if ( crange_min_glob(2) .lt. zero ) then ltheta = zero; utheta = pi lphi = zero; uphi = pi sym_factor = two center(2) = zero s_symy = .true. else if ( crange_min_glob(1) .lt. zero ) then ltheta = zero; utheta = pi lphi = -half * pi; uphi = half * pi sym_factor = two center(1) = zero s_symx = .true. end if end if ! Set the theta index to be used for the equatorial circumference integration ! based on the symmetries. if ( .not. s_symz ) then githeta = ( ntheta - 1 ) / 2 + 1 else githeta = ntheta end if ! Set the phi index to be used for the polar circumference integration based ! on the symmetries. if ( ( .not. s_symy ) .and. ( s_symx ) ) then gjphi = ( nphi - 1 ) / 2 + 1 else gjphi = 1 end if ! Set up the symmetry muliplication factors based on the range in the ! angular coordinates. theta_sym_factor = two * pi / ( utheta - ltheta ) phi_sym_factor = two * pi / ( uphi - lphi ) ! 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. dtheta = ( utheta - ltheta ) / ( ntheta - 1 ) dphi = ( uphi - lphi ) / ( nphi - 1 ) dthetainv = one / dtheta dphiinv = one / dphi do i = 1, lsh(1) ctheta(i,:) = ltheta + dtheta * ( i + lbnd(1) - 1 ) end do do j = 1, lsh(2) cphi(:,j) = lphi + dphi * ( j +lbnd(2) - 1 ) end do ! Calculate sines and cosines and store them in grid arrays since they ! are computationally expensive and memory for 2D grid arrays is cheap. sintheta = sin(ctheta) costheta = cos(ctheta) sinphi = sin(cphi) cosphi = cos(cphi) ! 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 ! and is initialised to 2*Delta. The grid array n_since_last_reduction is ! a help in keeping track on how long ago drsurf for a given point has been ! reduced and is initialised to -1 to indicate that no reduction has taken ! place. drsurf = two * min_delta rsurf = zero n_since_last_reduction = -1 ! Get and interpolation handle for required polynomial interpolation. ! If Lagrange polynomial interpolation does not work try Hermite ! polynomial interpolation, which avoids problems with non-continuous ! derivatives in the Newton iteration later on. ! First convert the string parameter to a Fortran string. call CCTK_FortranString ( surface_interp_len, surface_interpolator, & surface_interp ) ! Then get the handle ... call CCTK_InterpHandle ( interp_handle, surface_interp(1:surface_interp_len) ) ! The following is done because some C-preprocessors seem to have problems ! handling character strings that are continued on several lines. if ( interp_handle .lt. 0 ) then warn_message = 'Cannot get handle for interpolation. ' warn_message = trim(warn_message)//'Forgot to activate an implementation ' warn_message = trim(warn_message)//'providing interpolation operators?' call CCTK_WARN( 0, trim(warn_message) ) end if ! Write the required interpolation order parameter into a fortran string. write(surface_order,'(a6,i1)') 'order=',surface_interpolation_order ! Then create a table from the string. call Util_TableCreateFromString ( table_handle, surface_order ) if ( table_handle .lt. 0 ) then call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' ) end if ! Get the 3D coordinate system handle. call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' ) if ( coord_system_handle .lt. 0) then warn_message = 'Cannot get handle for cart3d coordinate system. ' warn_message = trim(warn_message)//'Forgot to activate an implementation ' warn_message = trim(warn_message)//'providing coordinates?' call CCTK_WARN( 0, trim(warn_message) ) endif ! Get pointers to the grid arrays containing the interpolation points. interp_coords(1) = CCTK_PointerTo(interp_x) interp_coords(2) = CCTK_PointerTo(interp_y) interp_coords(3) = CCTK_PointerTo(interp_z) ! Get pointer to the interpolation output array out_array(1) = CCTK_PointerTo(f_interp) ! Get the variable index for 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 ! The algorithm simply starts at the center and works it way outwards ! in all directions. When the function changes signs, it reduces the ! changes in radius until all points have negative f values and are ! within one grid cells difference from the surface. find_starting_points: do do j = 1, lsh(2) do i = 1, lsh(1) ! The radius that are tested in this direction. rloc = rsurf(i,j) + drsurf(i,j) ! Calculate the interpolation points. interp_x(i,j) = center(1) + rloc * sintheta(i,j) * cosphi(i,j) interp_y(i,j) = center(2) + rloc * sintheta(i,j) * sinphi(i,j) interp_z(i,j) = center(3) + rloc * costheta(i,j) ! If the point has been reduced once, increment n_since_last_reduction. 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 the interpolator. 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) ) 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 ! Find the maximal change in rsurf. maxdr_loc = maxval ( drsurf ) call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, & maxdr_loc, maxdr, CCTK_VARIABLE_REAL ) ! Find the maximal value of interpolated f values. maxf_loc = maxval ( abs ( f_interp ) ) call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, & maxf_loc, maxf, CCTK_VARIABLE_REAL ) ! If the maximal change in rsurf is small enough, we are satisfied that ! the Newton iteration will be able to find the root. if ( maxdr .lt. min_delta ) exit find_starting_points ! Check if we have crossed the zero point and if we have reduce the ! change in rsurf we are going to try next. 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 ! Now for the Newton iteration we also need the interpolates derivatives. 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 k = 1, 15 ! Calculate the interpolation points. do j = 1, lsh(2) do i = 1, lsh(1) interp_x(i,j) = center(1) + rsurf(i,j) * sintheta(i,j) * cosphi(i,j) interp_y(i,j) = center(2) + rsurf(i,j) * sintheta(i,j) * sinphi(i,j) interp_z(i,j) = center(3) + rsurf(i,j) * costheta(i,j) end do end do ! Call the interpolator. 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 ) 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 ! Find the maximum value of interpolated f values. maxf_loc = maxval ( abs ( f_interp ) ) call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, & maxf_loc, maxf, CCTK_VARIABLE_REAL ) ! If the maximum interpolated value of f is to big, perform a Newton ! iteration 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) * sintheta(i,j) * cosphi(i,j) + & dfdy_interp(i,j) * sintheta(i,j) * sinphi(i,j) + & dfdz_interp(i,j) * costheta(i,j) ) end do end do else ! Else we are done exit find_surface_points end if end do find_surface_points ! If we used to many iterations, the Newton iteration did not converge, so ! print a warning. if ( k .gt. 15 ) then warn_message = 'Newton iteration did not converge! Area may be inaccurate' call CCTK_WARN(1, trim(warn_message) ) end if drdtheta = zero drdphi = zero ! Calculate the derivatives of r. ! First for interior points do j = 2, lsh(2) - 1 do i = 2, lsh(1) - 1 drdtheta(i,j) = half * dthetainv * ( rsurf(i+1,j) - rsurf(i-1,j) ) drdphi(i,j) = half * dphiinv * ( rsurf(i,j+1) - rsurf(i,j-1) ) end do end do ! If boundary 1 is physical update those points. if ( bbox(1) .eq. 1 ) then do j = 2, lsh(2) - 1 drdtheta(1,j) = half * dthetainv * ( - three * rsurf(1,j) + & four * rsurf(2,j) - & rsurf(3,j) ) drdphi(1,j) = half * dphiinv * ( rsurf(1,j+1) - rsurf(1,j-1) ) end do end if ! If boundary 2 is physical update those points. if ( bbox(2) .eq. 1 ) then im = lsh(1) do j = 2, lsh(2) - 1 drdtheta(im,j) = half * dthetainv * ( three * rsurf(im,j) - & four * rsurf(im-1,j) + & rsurf(im-2,j) ) drdphi(im,j) = half * dphiinv * ( rsurf(im,j+1) - rsurf(im,j-1) ) end do end if ! If boundary 3 is physical update those points. if ( bbox(3) .eq. 1 ) then do i = 2, lsh(1) - 1 drdtheta(i,1) = half * dthetainv * ( rsurf(i+1,1) - rsurf(i-1,1) ) drdphi(i,1) = half * dphiinv * ( - three * rsurf(i,1) + & four * rsurf(i,2) - & rsurf(i,3) ) end do end if ! If boundary 4 is physical update those points. if ( bbox(4) .eq. 1 ) then jm = lsh(2) do i = 2, lsh(1) - 1 drdtheta(i,jm) = half * dthetainv * ( rsurf(i+1,jm) - rsurf(i-1,jm) ) drdphi(i,jm) = half * dphiinv * ( three * rsurf(i,jm) - & four * rsurf(i,jm-1) + & rsurf(i,jm-2) ) end do end if ! If corner 1 is physical update that point. if ( bbox(1) .eq. 1 .and. bbox(3) .eq. 1 ) then drdtheta(1,1) = half * dthetainv * ( - three * rsurf(1,1) + & four * rsurf(2,1) - & rsurf(3,1) ) drdphi(1,1) = half * dphiinv * ( - three * rsurf(1,1) + & four * rsurf(1,2) - & rsurf(1,3) ) end if ! If corner 2 is physical update that point. if ( bbox(1) .eq. 1 .and. bbox(4) .eq. 1 ) then jm = lsh(2) drdtheta(1,jm) = half * dthetainv * ( - three * rsurf(1,jm) + & four * rsurf(2,jm) - & rsurf(3,jm) ) drdphi(1,jm) = half * dphiinv * ( three * rsurf(1,jm) - & four * rsurf(1,jm-1) + & rsurf(1,jm-2) ) end if ! If corner 3 is physical update that point. if ( bbox(2) .eq. 1 .and. bbox(3) .eq. 1 ) then im = lsh(1) drdtheta(im,1) = half * dthetainv * ( three * rsurf(im,1) - & four * rsurf(im-1,1) + & rsurf(im-2,1) ) drdphi(im,1) = half * dphiinv * ( - three * rsurf(im,1) + & four * rsurf(im,2) - & rsurf(im,3) ) end if ! If corner 4 is physical update that point. if ( bbox(2) .eq. 1 .and. bbox(4) .eq. 1 ) then im = lsh(1) jm = lsh(2) drdtheta(im,jm) = half * dthetainv * ( three * rsurf(im,jm) - & four * rsurf(im-1,jm) + & rsurf(im-2,jm) ) drdphi(im,jm) = half * dphiinv * ( three * rsurf(im,jm) - & four * rsurf(im,jm-1) + & rsurf(im,jm-2) ) end if ! I do not remember why this is done. I will have to investigate. integrate_counter = integrate_counter - 1 ! Get rid of the interpolation table. call Util_TableDestroy ( status, table_handle ) 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 implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS more_surfaces = 1 surface_counter = 0 sc = zero end subroutine EHFinder_CountSurfacesInit subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS) use EHFinder_mod use EHFinder_fuzzy implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS CCTK_INT, dimension(3) :: minl CCTK_REAL :: minf_loc, minf_glob CCTK_INT :: procpos_loc, procpos_glob ! Find index ranges for points excluding ghost zones and symmetry points for ! the 3D grid functions. #include "include/physical_part.h" ! 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,levelset_counter), & fuzzy_eq ( sc(ixl:ixr,jyl:jyr,kzl:kzr), zero, small ) .and. & ( 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,levelset_counter) ! Find the minimum over all processors. call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & minf_loc, minf_glob, CCTK_VARIABLE_REAL ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of minf_glob failed') end if ! To figure out on which processor the minimum is located I set procpos_loc ! to the local processor number if the local minimum is equal to the global ! minimum otherwise I set it to the total number of processors. When I then ! find the minimum of procpos_loc over all processors I get the processor ! containing the minimum. if ( minf_loc .eq. minf_glob ) then procpos_loc = CCTK_MyProc ( cctkGH ) else procpos_loc = CCTK_nProcs ( cctkGH ) end if call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & procpos_loc, procpos_glob, CCTK_VARIABLE_INT ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of minf_glob failed') end if ! If minf_glob<0 there must exist another surface so lets set up the variables ! to start marking and counting these points. if ( minf_glob .lt. zero ) then surface_counter = surface_counter + 1 points_counter = 0 more_points = 1 if ( procpos_loc .eq. procpos_glob ) then sc(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1) = real(surface_counter,wp) end if else more_surfaces = 0 more_points = 0 end if end subroutine EHFinder_CountSurfaces ! This routine marks all points inside the current surface. subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) use EHFinder_mod use EHFinder_fuzzy implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS CCTK_INT :: n_loc, n_glob, i, j, k CCTK_INT :: i1, i2, j1, j2, k1, k2 logical :: marked ! Find index ranges for points excluding ghost zones and symmetry points for ! the 3D grid functions. #include "include/physical_part.h" ! Initialise the local point counter. n_loc = 0 ! Check if any points have a marked point as a neighbour. Only look at ! active points with f<0 that are not already marked. Take care at the ! boundaries (physical and excision) not to look at inactive or not ! exisiting points. do k = kzl, kzr do j = jyl, jyr do i = ixl, ixr marked = .false. if ( ( eh_mask(i,j,k,levelset_counter) .ge. 0 ) .and. & ( f(i,j,k,levelset_counter) .lt. 0 ) .and. & fuzzy_ne ( sc(i,j,k), real(surface_counter,wp), small ) ) then ! Find out if we are on the boundary if ( i + cctk_lbnd(1) .eq. 1 ) then i1 = 0 else i1 = -1 end if if ( i + cctk_lbnd(1) .eq. cctk_gsh(1) ) then i2 = 0 else i2 = 1 end if if ( j + cctk_lbnd(2) .eq. 1 ) then j1 = 0 else j1 = -1 end if if ( j + cctk_lbnd(2) .eq. cctk_gsh(2) ) then j2 = 0 else j2 = 1 end if if ( k + cctk_lbnd(3) .eq. 1 ) then k1 = 0 else k1 = -1 end if if ( k + cctk_lbnd(3) .eq. cctk_gsh(3) ) then k2 = 0 else k2 = 1 end if ! Then check if any valid neighbours are marked if ( any ( fuzzy_eq ( sc(i+i1:i+i2,j+j1:j+j2,k), & real(surface_counter,wp), small ) ) ) then marked = .true. end if if ( any ( fuzzy_eq ( sc(i+i1:i+i2,j,k+k1:k+k2), & real(surface_counter,wp), small ) ) ) then marked = .true. end if if ( any ( fuzzy_eq ( sc(i,j+j1:j+j2,k+k1:k+k2), & real(surface_counter,wp), small ) ) ) then marked = .true. end if if ( marked ) then sc(i,j,k) = real(surface_counter,wp) n_loc = n_loc + 1 end if end if end do end do end do ! Reduce the number of newly marked points over all processors call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & n_loc, n_glob, CCTK_VARIABLE_INT ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of n_glob failed') end if ! If the total number of newly marked points is 0, set more_points to 0 ! to indicate that we are finished. Otherwise continue and add the number ! of points to the points_counter. if ( n_glob .eq. 0 ) then more_points = 0 end if points_counter = points_counter + n_glob end subroutine EHFinder_MarkSurfaces ! This routine prints out information about the found surfaces. subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS) use EHFinder_mod implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS character(len=80) :: info_message ! Write out how many surfaces were found. if ( surface_counter .eq. 1 ) then write(info_message,'(a9,i3,a22,i3)') 'There is ',surface_counter, & ' surface in level set ',levelset_counter else write(info_message,'(a10,i3,a23,i3)') 'There are ',surface_counter, & ' surfaces in level set ',levelset_counter end if call CCTK_INFO( trim(info_message) ) integrate_counter = surface_counter end subroutine EHFinder_InfoSurfaces