diff options
| -rw-r--r-- | interface.ccl | 33 | ||||
| -rw-r--r-- | param.ccl | 3 | ||||
| -rw-r--r-- | schedule.ccl | 83 | ||||
| -rw-r--r-- | src/EHFinder_FindSurface.F90 | 40 | ||||
| -rw-r--r-- | src/EHFinder_Integrate2.F90 | 143 |
5 files changed, 232 insertions, 70 deletions
diff --git a/interface.ccl b/interface.ccl index f4fe515..9af0506 100644 --- a/interface.ccl +++ b/interface.ccl @@ -53,7 +53,8 @@ CCTK_INT surface_counter TYPE=SCALAR CCTK_INT points_counter TYPE=SCALAR CCTK_INT more_surfaces TYPE=SCALAR CCTK_INT more_points TYPE=SCALAR - +CCTK_INT integrate_counter TYPE=SCALAR + 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, sintheta, costheta, sinphi, cosphi, @@ -90,9 +91,31 @@ 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_REAL eh_area2 TYPE=SCALAR -CCTK_REAL eh_centroid TYPE=SCALAR + +#CCTK_REAL eh_area TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT + +CCTK_REAL eh_area TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT + +#CCTK_REAL eh_area[maximum_surface_number] TYPE=SCALAR + +#CCTK_REAL eh_area2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT #"The area of the event horizon(s)" +CCTK_REAL eh_area2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT #"The area of the event horizon(s)" + +#CCTK_REAL eh_area2[maximum_surface_number] TYPE=SCALAR +#CCTK_REAL eh_area3 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT + +CCTK_REAL eh_centroid TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT { eh_centroid_x, eh_centroid_y, eh_centroid_z -} "The centroid of the event horizon" +} "The centroid of the event horizon(s)" + +CCTK_REAL eh_centroid2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT +{ + eh_centroid2_x, eh_centroid2_y, eh_centroid2_z +} "The centroid of the event horizon(s)" + +#CCTK_REAL eh_centroid3 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT +#{ +# eh_centroid3_x, eh_centroid3_y, eh_centroid3_z +#} "The centroid of the event horizon(s)" + @@ -178,7 +178,8 @@ INT n_array_ghosts "Number of ghost points in the surface grid array" 1: :: "Positive please" } 1 -INT maximum_surface_number "The maximum number of surfaces expected in the data"{ +INT maximum_surface_number "The maximum number of surfaces expected in the data" +{ 1:* :: "Positive please" } 1 diff --git a/schedule.ccl b/schedule.ccl index 572ea93..7652251 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -17,7 +17,8 @@ else STORAGE: mask_bak STORAGE: surface_arrays STORAGE: center_arrays - STORAGE: eh_area2 +# STORAGE: eh_area, eh_centroid + STORAGE: eh_area2, eh_centroid2 # STORAGE: surface_index } @@ -88,77 +89,99 @@ schedule EHFinder_InitWeights at CCTK_POSTINITIAL schedule GROUP EHFinder_Surfaces at CCTK_ANALYSIS { STORAGE: more_surfaces, surface_counter, more_points, points_counter + STORAGE: integrate_counter STORAGE: surface_index - TRIGGER: eh_area2, eh_centroid -} "Schedule group for counting the number of surfaces and integrating over them" +# STORAGE: eh_area2, eh_centroid2 + TRIGGERS: eh_area, eh_centroid +} "Count the number of surfaces and integrate over them" schedule EHFinder_CountSurfacesInit in EHFinder_Surfaces { LANG: Fortran - TRIGGER: eh_area2, eh_centroid -} "Initialising while loop control for the group EHFinder_CountMarkSurfaces" + TRIGGERS: eh_area, eh_centroid +} "Initialise while loop control" schedule GROUP EHFinder_CountMarkSurfaces in EHFinder_Surfaces after EHFinder_CountSurfacesInit WHILE ehfinder::more_surfaces { - TRIGGER: eh_area2, eh_centroid -} "Scedule group for counting and marking surfaces" + TRIGGERS: eh_area, eh_centroid +} "Counting and mark surfaces" schedule EHFinder_CountSurfaces in EHFinder_CountMarkSurfaces { LANG: Fortran - TRIGGER: eh_area2, eh_centroid + TRIGGERS: eh_area, eh_centroid SYNC: surface_index -} "Check if there are more surfaces and initialise while loop control for EHFinder_MarkSurfaces" +} "Check if there are more surfaces" schedule EHFinder_MarkSurfaces in EHFinder_CountMarkSurfaces after EHFinder_CountSurfaces WHILE ehfinder::more_points { LANG: Fortran - TRIGGER: eh_area2, eh_centroid + TRIGGERS: eh_area, eh_centroid SYNC: surface_index -} "Mark points inside the current surface until they are all marked" +} "Mark points inside the current surface" schedule EHFinder_InfoSurfaces in EHFinder_Surfaces after EHFinder_CountMarkSurfaces { LANG: Fortran - TRIGGER: eh_area2, eh_centroid -} "Output input about found surfaces" + TRIGGERS: eh_area, eh_centroid +} "Output info about found surfaces" -schedule EHFinder_FindSurface at CCTK_ANALYSIS +schedule group EHFinder_Integration in EHFinder_Surfaces after EHFinder_InfoSurfaces while ehfinder::integrate_counter { - LANG: Fortran - STORAGE: surface_tmp_arrays, surface_int_array - TRIGGER: eh_area2, eh_centroid -} "Find Surface" + TRIGGERS: eh_area, eh_centroid +} "Find and integrate over surfaces" -schedule EHFinder_Integrate at CCTK_ANALYSIS after EHFinder_FindSurface +schedule EHFinder_FindSurface in EHFinder_Integration { LANG: Fortran - STORAGE: eh_area - TRIGGER: eh_area -} "Calculate integrals" + STORAGE: surface_tmp_arrays, surface_int_array + TRIGGERS: eh_area, eh_centroid +} "Find Surface" -schedule EHFinder_FindSurfaceElement at CCTK_ANALYSIS after EHFinder_FindSurface +schedule EHFinder_FindSurfaceElement in EHFinder_Integration after EHFinder_FindSurface { LANG: Fortran STORAGE: surface_tmp_arrays, interp_metric_arrays - TRIGGER: eh_area2, eh_centroid + TRIGGERS: eh_area, eh_centroid } "Find Surface Area Element" -schedule EHFinder_IntegrateArea at CCTK_ANALYSIS after EHFinder_FindSurfaceElement +schedule EHFinder_IntegrateArea in EHFinder_Integration after EHFinder_FindSurfaceElement { LANG: Fortran - STORAGE: integrate_tmp_array - TRIGGER: eh_area2, eh_centroid + STORAGE: integrate_tmp_array + TRIGGERS: eh_area, eh_centroid } "Calculate area integrals" -schedule EHFinder_IntegrateCentroid at CCTK_ANALYSIS after EHFinder_IntegrateArea +schedule EHFinder_IntegrateCentroid in EHFinder_Integration after EHFinder_IntegrateArea { LANG: Fortran - STORAGE: eh_centroid, surface_tmp_arrays, integrate_tmp_array - TRIGGER: eh_centroid + STORAGE: surface_tmp_arrays, integrate_tmp_array + TRIGGERS: eh_centroid } "Calculate centroid integrals" +schedule EHFinder_CopyArea at CCTK_ANALYSIS after EHFinder_Surfaces +{ + LANG: Fortran + STORAGE: eh_area + TRIGGERS: eh_area +} "Copy areas to output variable" + +schedule EHFinder_CopyCentroid at CCTK_ANALYSIS after EHFinder_Surfaces +{ + LANG: Fortran + STORAGE: eh_centroid + TRIGGERS: eh_centroid +} "Copy areas to output variable" + +#schedule EHFinder_Integrate at CCTK_ANALYSIS +#schedule EHFinder_Integrate at CCTK_ANALYSIS +#{ +# LANG: Fortran +# STORAGE: eh_area +# TRIGGERS: eh_area +#} "Calculate integrals" + # Read in the data used in reconstructing the 4-metric if necessary if (CCTK_Equals(eh_metric_type,"numerical")) diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90 index fe41fd1..41bbac4 100644 --- a/src/EHFinder_FindSurface.F90 +++ b/src/EHFinder_FindSurface.F90 @@ -15,11 +15,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k, im, jm + CCTK_INT :: i, j, k, 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=30) :: info_message CCTK_INT, dimension(4) :: bbox CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost, maxl @@ -42,7 +43,10 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) ! the 3D grid functions. #include "include/physical_part.h" - call CCTK_INFO ( 'Finding points on surface' ) + sn = surface_counter - integrate_counter + 1 + + write(info_message,'(a26,i3)') 'Finding points on surface ',sn + call CCTK_INFO ( info_message ) ! Find information about the local and global properties of the 2D ! Grid arrays. @@ -77,26 +81,33 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if ! If the user defined center point should be used... - if ( use_user_center ) then + 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. + ! are within one grid cell distance from the surface in question. 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 + 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 ) ) .and. & + ( ( sc(i,j,k) .eq. sn ) .or. & + ( sc(i-1,j,k) .eq. sn ) .or. & + ( sc(i+1,j,k) .eq. sn ) .or. & + ( sc(i,j-1,k) .eq. sn ) .or. & + ( sc(i,j+1,k) .eq. sn ) .or. & + ( sc(i,j,k-1) .eq. sn ) .or. & + ( sc(i,j,k+1) .eq. sn ) ) ) 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) @@ -115,7 +126,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) CCTK_VARIABLE_REAL ) ncinv = one / nc center = ncinv * center -! print*,center +! print*,'EHFinder_FindSurfaces : ', sn, center ! For now symmetries are not handled properly if ( symx ) center(1) = zero @@ -183,8 +194,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) ! 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 ??" ) + call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. Forgot to activate an implementation providing coordinates ??" ) endif ! Get pointers to the grid arrays containing the interpolation points. @@ -427,6 +437,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" ) + integrate_counter = integrate_counter - 1 ! print*,'Calculated derivatives' ! print*,rsurf call Util_TableDestroy ( status, table_handle ) @@ -634,4 +645,7 @@ subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS) end if call CCTK_INFO(info_message) + integrate_counter = surface_counter + +! print*,'EHFinder_InfoSurfaces : ',integrate_counter end subroutine EHFinder_InfoSurfaces diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90 index f82cd54..e463a43 100644 --- a/src/EHFinder_Integrate2.F90 +++ b/src/EHFinder_Integrate2.F90 @@ -72,8 +72,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) 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??" ) + 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" ) @@ -83,8 +82,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) 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 ??" ) + 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) @@ -201,10 +199,14 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: int_var + CCTK_INT :: int_var, sn + CCTK_REAL :: eh_area_tmp character(len=30) :: info_message +! CCTK_INT, dimension(1) :: lbnd, ubnd, lsh call CCTK_INFO ( 'Integrating area' ) + + sn = surface_counter - integrate_counter int_tmp = weights * da call CCTK_VarIndex ( int_var, "ehfinder::int_tmp" ) @@ -212,10 +214,33 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) Call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' ) end if call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & - eh_area2, 1, int_var ) - write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2 + eh_area2(sn), 1, int_var ) + write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2(sn) + +! call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & +! eh_area_tmp, 1, int_var ) +! write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area_tmp + call CCTK_INFO(info_message) +! call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_area2" ) +! if ( ierr .lt. 0 ) then +! call CCTK_WARN ( 0, "cannot get lower bounds for area array" ) +! end if +! call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_area2" ) +! if ( ierr .lt. 0 ) then +! call CCTK_WARN ( 0, "cannot get upper bounds for area array" ) +! end if +! call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_area2" ) +! if ( ierr .lt. 0 ) then +! call CCTK_WARN ( 0, "cannot get local size for area array" ) +! end if +! +! if ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then +! eh_area2(sn-lbnd(1)) = eh_area_tmp +! end if +! print*,lbnd,ubnd,lsh + end subroutine EHFinder_IntegrateArea @@ -229,10 +254,15 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: int_var + CCTK_INT :: int_var, sn, k character(len=50) :: info_message +! CCTK_INT, dimension(1) :: lbnd, ubnd, lsh + CCTK_REAL :: centroid_x, centroid_y, centroid_z call CCTK_INFO ( 'Integrating centroid' ) + + sn = surface_counter - integrate_counter + interp_x = center(1) + rsurf * sintheta * cosphi interp_y = center(2) + rsurf * sintheta * sinphi interp_z = center(3) + rsurf * costheta @@ -244,31 +274,66 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) int_tmp = interp_x * weights * da call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & - eh_centroid_x, 1, int_var ) + eh_centroid2_x(sn), 1, int_var ) int_tmp = interp_y * weights * da call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & - eh_centroid_y, 1, int_var ) + eh_centroid2_y(sn), 1, int_var ) int_tmp = interp_z * weights * da call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & - eh_centroid_z, 1, int_var ) - - eh_centroid_x = eh_centroid_x / eh_area2 - eh_centroid_y = eh_centroid_y / eh_area2 - eh_centroid_z = eh_centroid_z / eh_area2 - - write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ',eh_centroid_x, & - eh_centroid_y, & - eh_centroid_z + eh_centroid2_z(sn), 1, int_var ) + +! call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & +! centroid_x, 1, int_var ) +! +! int_tmp = interp_y * weights * da +! call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & +! centroid_y, 1, int_var ) +! +! int_tmp = interp_z * weights * da +! call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & +! centroid_z, 1, int_var ) +! +! call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_area2" ) +! if ( ierr .lt. 0 ) then +! call CCTK_WARN ( 0, "cannot get lower bounds for area array" ) +! end if +! call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_area2" ) +! if ( ierr .lt. 0 ) then +! call CCTK_WARN ( 0, "cannot get upper bounds for area array" ) +! end if +! call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_area2" ) +! if ( ierr .lt. 0 ) then +! call CCTK_WARN ( 0, "cannot get local size for area array" ) +! end if + + eh_centroid2_x(sn) = eh_centroid2_x(sn) / eh_area2(sn) + eh_centroid2_y(sn) = eh_centroid2_y(sn) / eh_area2(sn) + eh_centroid2_z(sn) = eh_centroid2_z(sn) / eh_area2(sn) + +! if ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then +! k = sn - lbnd(1) +! eh_centroid2_x(k) = eh_centroid2_x(k) / eh_area2(k) +! eh_centroid2_y(k) = eh_centroid2_y(k) / eh_area2(k) +! eh_centroid2_z(k) = eh_centroid2_z(k) / eh_area2(k) +! end if + + write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', & + eh_centroid2_x(sn), & + eh_centroid2_y(sn), & + eh_centroid2_z(sn) + +! write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', & +! centroid_x, & +! centroid_y, & +! centroid_z call CCTK_INFO(info_message) end subroutine EHFinder_IntegrateCentroid ! This routine sets up the weights for the Simpsons rule integration - -! This routine sets up the weights for the Simpsons rule integration ! over the surface. subroutine EHFinder_InitWeights(CCTK_ARGUMENTS) @@ -340,3 +405,39 @@ subroutine EHFinder_InitWeights(CCTK_ARGUMENTS) ! 4/9 16/9 8/9 16/9 8/9 16/9 4/9 ! 1/9 4/9 2/9 4/9 2/9 4/9 1/9 end subroutine EHFinder_InitWeights + + +subroutine EHFinder_CopyArea(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + call CCTK_INFO('Copying area data') + + eh_area = eh_area2 + +end subroutine EHFinder_CopyArea + + +subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + 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 subroutine EHFinder_CopyCentroid |