diff options
Diffstat (limited to 'src/EHFinder_FindSurface.F90')
| -rw-r--r-- | src/EHFinder_FindSurface.F90 | 259 |
1 files changed, 136 insertions, 123 deletions
diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90 index 2537408..d53db6e 100644 --- a/src/EHFinder_FindSurface.F90 +++ b/src/EHFinder_FindSurface.F90 @@ -21,13 +21,14 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) 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, maxl + CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost CCTK_REAL, dimension(3) :: center_loc CCTK_INT :: nc_loc, nc @@ -55,52 +56,52 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) ! 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 ( 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" ) +! 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" ) + 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" ) + 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" ) + 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" ) + 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" ) + 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" ) + 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" ) + 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" ) + 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" ) + 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" ) + 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" ) + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) end if - ! If the user defined center point should be used... +! 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. +! 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 @@ -130,7 +131,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end do end do - ! Reduce over the processors. +! 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, & @@ -138,14 +139,10 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) CCTK_VARIABLE_REAL ) ncinv = one / nc center = ncinv * center -! print*,'EHFinder_FindSurfaces : ', sn, center -! print* end if -! print*,'Center : ',center -! print*,'use_user_center : ',use_user_center - ! Try to figure out the symmetries. First find the range in coordinates - ! for the points inside the current surface. +! 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 = sc .eq. sn ) crange_min_loc(2) = minval ( y, mask = sc .eq. sn ) crange_min_loc(3) = minval ( z, mask = sc .eq. sn ) @@ -153,18 +150,15 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) crange_max_loc(2) = maxval ( y, mask = sc .eq. sn ) crange_max_loc(3) = maxval ( z, mask = sc .eq. sn ) +! 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 ) -! print*,crange_min_loc,crange_max_loc -! print*,crange_min_glob,crange_max_glob -! print* -! print*,'Symmetries : ', symx, symy, symz - ! Full mode. Modify if necessary. +! Full mode. Modify if necessary. ltheta = zero; utheta = pi lphi = zero; uphi = two * pi sym_factor = one @@ -172,9 +166,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) if ( CCTK_EQUALS ( domain, 'bitant' ) ) then - ! Bitant with symmetry in the z-direction +! 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 @@ -184,7 +181,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Bitant with symmetry in the y-direction +! Bitant with symmetry in the y-direction if ( CCTK_EQUALS ( bitant_plane, 'xz' ) ) then if ( crange_min_glob(2) .lt. zero ) then @@ -196,7 +193,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Bitant with symmetry in the x-direction +! Bitant with symmetry in the x-direction if ( CCTK_EQUALS ( bitant_plane, 'yz' ) ) then if ( crange_min_glob(1) .lt. zero ) then @@ -211,7 +208,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) if ( CCTK_EQUALS ( domain, 'quadrant' ) ) then - ! Quadrant in x-direction +! Quadrant in x-direction if ( CCTK_EQUALS ( quadrant_direction, 'x' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & @@ -236,7 +233,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Quadrant in y-direction +! Quadrant in y-direction if ( CCTK_EQUALS ( quadrant_direction, 'y' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & @@ -261,7 +258,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Quadrant in z-direction +! Quadrant in z-direction if ( CCTK_EQUALS ( quadrant_direction,'z' ) ) then if ( ( crange_min_glob(2) .lt. zero ) .and. & @@ -287,7 +284,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Octant +! Octant if ( CCTK_EQUALS ( domain, 'octant' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & @@ -340,24 +337,30 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) 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. +! 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 ) @@ -370,64 +373,71 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) cphi(:,j) = lphi + dphi * ( j +lbnd(2) - 1 ) end do -! print*,'Delta : ',dtheta,dphi - ! Calculate sines and cosines and store them in grid arrays since they - ! are expensive. +! 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. +! 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. +! 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. +! First convert the string parameter to a Fortran string. call CCTK_FortranString ( surface_interp_len, surface_interpolator, & surface_interp ) - ! Then get the hande ... +! 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 - call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" ) + warn_message = 'Cannot get handle for interpolation. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing interpolation operators??' + call CCTK_WARN( 0, warn_message ) end if - ! Write the required interpolation order parameter into a fortran string. +! 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. +! 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" ) + 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" ) +! Get the 3D coordinate system handle. + 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 ??" ) + warn_message = 'Cannot get handle for cart3d coordinate system. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing coordinates??' + call CCTK_WARN( 0, warn_message ) endif - ! Get pointers to the grid arrays containing the interpolation points. +! 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 +! Get pointer to the interpolation output array out_array(1) = CCTK_PointerTo(f_interp) - ! Get the variable index for f +! Get the variable index for f write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1, ']' call CCTK_VarIndex ( in_array, vname ) @@ -442,22 +452,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) do j = 1, lsh(2) do i = 1, lsh(1) - ! The radius that are tested in this direction. +! The radius that are tested in this direction. rloc = rsurf(i,j) + drsurf(i,j) - ! Calculate the interpolation points. +! 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 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, & @@ -470,18 +480,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) 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 ) - maxl = maxloc ( abs ( f_interp ) ) - +! 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 @@ -497,28 +511,28 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end do end do find_starting_points -! print*,'Max values: ',maxval(rsurf),maxval(f_interp) -! print*,'Min values: ',minval(rsurf),minval(f_interp) -! print*,'found 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" ) + op_indices, 'operand_indices' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + call CCTK_WARN ( 0, 'Cannot set operand indices array in parameter table' ) endif call Util_TableSetIntArray ( status, table_handle, 4, & - op_codes, "operation_codes" ) + op_codes, 'operation_codes' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + 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) @@ -527,6 +541,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) 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, & @@ -539,12 +554,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) 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 ) - maxl = maxloc ( abs ( f_interp ) ) - +! 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) @@ -555,16 +571,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end do end do else +! Else we are done exit find_surface_points end if -! call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" ) end do find_surface_points - if ( k .gt. 15 ) call CCTK_WARN(1,'Newton iteration did not converge! Area may be inaccurate') -! print*,'found_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, warn_message ) + end if drdtheta = zero drdphi = zero + ! Calculate the derivatives of r. ! First for interior points @@ -661,13 +683,11 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) rsurf(im,jm-2) ) end if - call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" ) - +! I do not remember why this is done. I will have to investigate. integrate_counter = integrate_counter - 1 -! print*,'Calculated derivatives' -! print*,rsurf + +! Get rid of the interpolation table. call Util_TableDestroy ( status, table_handle ) -! print*,'gjphi = ',gjphi end subroutine EHFinder_FindSurface @@ -718,7 +738,6 @@ subroutine EHFinder_CountSurfacesInit(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS -! print*,'EHFinder_CountSurfacesInit called' more_surfaces = 1 surface_counter = 0 sc = 0 @@ -744,27 +763,27 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS) ! 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. +! 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), & ( sc(ixl:ixr,jyl:jyr,kzl:kzr) .eq. 0 ) .and. & ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) ) - ! Find the value of f at that location. +! 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. +! 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. +! 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 @@ -777,8 +796,8 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS) 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<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 @@ -813,13 +832,13 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) ! the 3D grid functions. #include "include/physical_part.h" - ! Initialise the local point counter. +! 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. +! 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 @@ -828,7 +847,7 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) ( 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 +! Find out if we are on the boundary if ( i + cctk_lbnd(1) .eq. 1 ) then i1 = 0 else @@ -881,31 +900,25 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) end do end do - ! Reduce the number of newly marked points over all processors +! 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 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 - ! This should not be necessary, since I have a SYNC in the schedule.ccl, - ! but for some reason that does not work. - call CCTK_SyncGroup ( ierr, cctkGH, "ehfinder::surface_index" ) - - call CartSymVN(ierr,cctkGH,'ehfinder::sc') - end subroutine EHFinder_MarkSurfaces -!This routine prints out information about the found surfaces. +! This routine prints out information about the found surfaces. subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS) use EHFinder_mod @@ -918,7 +931,7 @@ subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS) character(len=80) :: info_message - ! Write out how many surfaces were found. +! 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 |