diff options
Diffstat (limited to 'src/EHFinder_Integrate.F90')
-rw-r--r-- | src/EHFinder_Integrate.F90 | 704 |
1 files changed, 704 insertions, 0 deletions
diff --git a/src/EHFinder_Integrate.F90 b/src/EHFinder_Integrate.F90 new file mode 100644 index 0000000..db2a315 --- /dev/null +++ b/src/EHFinder_Integrate.F90 @@ -0,0 +1,704 @@ +! Integrating over the surface(s). Right now only in full mode. +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, l, im, jm + CCTK_INT :: interp_handle, table_handle, coord_system_handle + + character(len=200) :: area_interp + character(len=128) :: warn_message + CCTK_INT :: area_interp_len + character(len=7) :: area_order + + CCTK_INT, dimension(4) :: bbox + CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost + + CCTK_REAL :: dtheta, dphi, dthetainv, dphiinv + CCTK_REAL :: dxdth, dxdph, dydth, dydph, dzdth, dzdph + CCTK_POINTER, dimension(3) :: interp_coords + CCTK_POINTER, dimension(7) :: out_array + CCTK_INT, dimension(7) :: in_array + CCTK_INT, dimension(7), parameter :: out_types = (/ CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL /) + +! If finding of surface failed do not try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + +! Store the levelset_counter in a shorter named variable for convenience. + l = levelset_counter + +! Write an Info message about what we are doing. + call CCTK_INFO ( 'Finding surface element' ) + +! Obtain the dimensions and characteristics of the grid arrays. + call CCTK_GroupbboxGN ( ierr, cctkGH, 4, bbox, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get bounding box for surface arrays' ) + end if + call CCTK_GroupgshGN ( ierr, cctkGH, 2, gsh, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get global size for surface arrays' ) + end if + call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get lower bounds for surface arrays' ) + end if + call CCTK_GroupubndGN ( ierr, cctkGH, 2, ubnd, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get upper bounds for surface arrays' ) + end if + call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) + end if + call CCTK_GroupnghostzonesGN ( ierr, cctkGH, 2, nghost, & + 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) + end if + +! Setup the theta and phi spacings and their inverse. + dtheta = ctheta(2,1) - ctheta(1,1) + dphi = cphi(1,2) - cphi(1,1) + dthetainv = one / dtheta + dphiinv = one / dphi + +! Convert the interpolator parameter into a fortran string. + call CCTK_FortranString ( area_interp_len, area_interpolator, & + area_interp ) + +! Get an interpolation handle. + call CCTK_InterpHandle ( interp_handle, area_interp(1:area_interp_len) ) + + if ( interp_handle .lt. 0 ) then + 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 interpolation order parameter into a string + write(area_order,'(a6,i1)') 'order=',area_interpolation_order + +! Create a table from the interpolation order string. + call Util_TableCreateFromString ( table_handle, area_order ) + if ( table_handle .lt. 0 ) then + call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' ) + end if + +! Get a handle for the coordinate system. + 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 = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing coordinates??' + call CCTK_WARN( 0, warn_message ) + endif + +! Get the pointers to 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 the pointers to the interpolation return arrays. + out_array(1) = CCTK_PointerTo(gxxi) + out_array(2) = CCTK_PointerTo(gxyi) + out_array(3) = CCTK_PointerTo(gxzi) + out_array(4) = CCTK_PointerTo(gyyi) + out_array(5) = CCTK_PointerTo(gyzi) + out_array(6) = CCTK_PointerTo(gzzi) + out_array(7) = CCTK_PointerTo(psii) + +! find the cartesian coordinates for 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 + +! Get the indices for the grid functions to be interpolated. + call CCTK_VarIndex ( in_array(1), 'admbase::gxx' ) + call CCTK_VarIndex ( in_array(2), 'admbase::gxy' ) + call CCTK_VarIndex ( in_array(3), 'admbase::gxz' ) + call CCTK_VarIndex ( in_array(4), 'admbase::gyy' ) + call CCTK_VarIndex ( in_array(5), 'admbase::gyz' ) + call CCTK_VarIndex ( in_array(6), 'admbase::gzz' ) + +! If the static conformal factor is used, we also need to interpolate it. + if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then + call CCTK_VarIndex ( in_array(7), 'staticconformal::psi' ) + +! Call the interpolator for the metric and the conformal factor. + call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1) * lsh(2), CCTK_VARIABLE_REAL, & + interp_coords, 7, in_array, & + 7, out_types, out_array ) + +! Get the physical metric. + gxxi = psii**4 * gxxi + gxyi = psii**4 * gxyi + gxzi = psii**4 * gxzi + gyyi = psii**4 * gyyi + gyzi = psii**4 * gyzi + gzzi = psii**4 * gzzi + else +! Else call the interpolater only for the metric. + call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1) * lsh(2), CCTK_VARIABLE_REAL, & + interp_coords, 6, in_array(1:6), & + 6, out_types(1:6), out_array(1:6) ) + end if + +! For all points on the surface... + do j = 1, lsh(2) + do i = 1, lsh(1) + +! Calculate the coordinate transformation between the 3-d cartesian +! coordinates and the 2-d (theta,phi) coordinates on the surface. + dxdth = ( drdtheta(i,j) * sintheta(i,j) + & + rsurf(i,j) * costheta(i,j) ) * cosphi(i,j) + dydth = ( drdtheta(i,j) * sintheta(i,j) + & + rsurf(i,j) * costheta(i,j) ) * sinphi(i,j) + dzdth = drdtheta(i,j) * costheta(i,j) - rsurf(i,j) * sintheta(i,j) + dxdph = ( drdphi(i,j) * cosphi(i,j) - & + rsurf(i,j) * sinphi(i,j) ) * sintheta(i,j) + dydph = ( drdphi(i,j) * sinphi(i,j) + & + rsurf(i,j) * cosphi(i,j) ) * sintheta(i,j) + dzdph = drdphi(i,j) * costheta(i,j) + +! Using this coordinate transformation, calculate the induced 2-metric +! on the surface. + gtt(i,j) = dxdth**2 * gxxi(i,j) + dydth**2 * gyyi(i,j) + & + dzdth**2 * gzzi(i,j) + & + two * ( dxdth * dydth * gxyi(i,j) + & + dxdth * dzdth * gxzi(i,j) + & + dydth * dzdth * gyzi(i,j) ) + gtp(i,j) = dxdth * dxdph * gxxi(i,j) + & + ( dxdth * dydph + dydth * dxdph ) * gxyi(i,j) + & + ( dxdth * dzdph + dzdth * dxdph ) * gxzi(i,j) + & + dydth * dydph * gyyi(i,j) + & + ( dydth * dzdph + dzdth * dydph ) * gyzi(i,j) + & + dzdth * dzdph * gzzi(i,j) + gpp(i,j) = dxdph**2 * gxxi(i,j) + dydph**2 * gyyi(i,j) + & + dzdph**2 * gzzi(i,j) + & + two * ( dxdph * dydph * gxyi(i,j) + & + dxdph * dzdph * gxzi(i,j) + & + dydph * dzdph * gyzi(i,j) ) + +! Calculate the area element as the square-root of the determinant +! of the two metric multiplied by dtheta*dphi. + da(i,j) = dtheta * dphi * sqrt ( gtt(i,j) * gpp(i,j) - gtp(i,j)**2 ) + +! Calculate the arc length element for the polar circumference integration. + dltheta(i,j) = dtheta * sqrt ( gtt(i,j) ) + +! Calculate the arc length element for the equatorial circumference +! integration. + dlphi(i,j) = dphi * sqrt ( gpp(i,j) ) + + end do + end do + +end subroutine EHFinder_FindSurfaceElement + + +subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: int_var, sn, l + CCTK_REAL :: eh_area_tmp + character(len=30) :: info_message + CCTK_INT, dimension(1) :: lbnd, ubnd, lsh + +! If finding of surface failed do not try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + +! Store the levelset_counter in a shorter named variable for convenience. + l = levelset_counter + +! Write an info message, to indicate what we are doing. + call CCTK_INFO ( 'Integrating area' ) + + sn = surface_counter - integrate_counter + +! Taking into account the sym_factor and the weights, the area is simply +! The sum of the elements of this array. + int_tmp = sym_factor * weights * da + +! Get the index to the area integration array. + call CCTK_VarIndex ( int_var, 'ehfinder::int_tmp' ) + if ( int_var .lt. 0 ) then + call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' ) + end if + +! Do a sum reduce over all processors. The result is the area. + call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & + eh_area_tmp, 1, int_var ) + +! Write an info message reporting the area found. + write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area_tmp + + call CCTK_INFO(info_message) + +! Figure out the bounds of the area storage grid array. + 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 this processor contains the right part of the area storage grid array +! save the area, otherwise do nothing. + if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then + eh_area2(sn-lbnd(1),l) = eh_area_tmp + end if + +end subroutine EHFinder_IntegrateArea + + +subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: int_var, sn, k, l + character(len=50) :: info_message + CCTK_INT, dimension(1) :: lbnd, ubnd, lsh + CCTK_REAL :: centroid_x, centroid_y, centroid_z + +! If finding of surface failed do not try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + +! Store the levelset_counter in a shorter named variable for convenience. + l = levelset_counter + +! Write an info message, to indicate what we are doing. + call CCTK_INFO ( 'Integrating centroid' ) + + sn = surface_counter - integrate_counter + +! Calculate the (x,y,z) coordinates of the points on the surface. + interp_x = center(1) + rsurf * sintheta * cosphi + interp_y = center(2) + rsurf * sintheta * sinphi + interp_z = center(3) + rsurf * costheta + +! Get the index to the temporary integration array. + call CCTK_VarIndex ( int_var, 'ehfinder::int_tmp' ) + if ( int_var .lt. 0 ) then + call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' ) + end if + +! Take into account the symmetry, weights, x-coordinate and area element +! to calculate the integral of x*da over the surface. This gives +! centroid_x * Area, so later we divide by the area. + int_tmp = sym_factor * interp_x * weights * da + +! Sum up the all elements of int_tmp. + call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & + centroid_x, 1, int_var ) + +! Take into account the symmetry, weights, y-coordinate and area element +! to calculate the integral of y*da over the surface. This gives +! centroid_y * Area, so later we divide by the area. + int_tmp = sym_factor * interp_y * weights * da + +! Sum up the all elements of int_tmp. + call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & + centroid_y, 1, int_var ) + +! Take into account the symmetry, weights, z-coordinate and area element +! to calculate the integral of z*da over the surface. This gives +! centroid_z * Area, so later we divide by the area. + int_tmp = sym_factor * interp_z * weights * da + +! Sum up the all elements of int_tmp. + call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & + centroid_z, 1, int_var ) +! If this surface has symmetries set the corresponding centroids to zero. + if ( s_symx ) centroid_x = zero + if ( s_symy ) centroid_y = zero + if ( s_symz ) centroid_z = zero + +! Get the bounds for the centroid arrays. It should be enough to do +! it just for the x-array, since they are all defined in the same way. + call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, 'ehfinder::eh_centroid2_x' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get lower bounds for centroid array' ) + end if + call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, 'ehfinder::eh_centroid2_x' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get upper bounds for centroid array' ) + end if + call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, 'ehfinder::eh_centroid2_x' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for centroid array' ) + end if + +! If we are on the right processor, store the centroid, otherwise do nothing. + if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then + k = sn - lbnd(1) + eh_centroid2_x(k,l) = centroid_x / eh_area2(k,l) + eh_centroid2_y(k,l) = centroid_y / eh_area2(k,l) + eh_centroid2_z(k,l) = centroid_z / eh_area2(k,l) + end if + +! Write an info message with the calculated centroids. + write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', & + eh_centroid2_x(sn,l), & + eh_centroid2_y(sn,l), & + eh_centroid2_z(sn,l) + call CCTK_INFO(info_message) + +end subroutine EHFinder_IntegrateCentroid + + +subroutine EHFinder_IntegrateCircumference(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: int_var, sn, k, l + character(len=50) :: info_message + CCTK_INT, dimension(1) :: lbnd1, ubnd1, lsh1 + CCTK_INT, dimension(2) :: lbnd, lsh, nghost + CCTK_INT, dimension(4) :: bbox + 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 do not try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + +! Store the levelset_counter in a shorter named variable for convenience. + l = levelset_counter + + sn = surface_counter - integrate_counter + + ! Find out the lower bounds of the distributed integration grid arrays. + call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get lower bounds for surface arrays' ) + end if + + ! Find out the local size of the distributed integration grid arrays + call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) + end if + + ! Find out the bounding box of the distributed integration grid arrays + call CCTK_GroupbboxGN ( ierr, cctkGH, 4, bbox, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get bounding box for surface arrays' ) + end if + + ! Find out the ghost size of the distributed integration grid arrays + call CCTK_GroupnghostzonesGN ( ierr, cctkGH, 2, nghost, & + 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get ghost zones for surface arrays' ) + end if + +! Initialize the lower and upper local index for the theta and phi arrays. + ithl = 1; ithr = lsh(1) + jphl = 1; jphr = lsh(2) + +! If the boundaries are internal, adjust the lower and upper local index +! to take into account the ghost cells. + if ( bbox(1) .eq. 0 ) ithl = ithl + nghost(1) + if ( bbox(2) .eq. 0 ) ithr = ithr - nghost(1) + if ( bbox(3) .eq. 0 ) jphl = jphl + nghost(2) + if ( bbox(4) .eq. 0 ) jphr = jphr - nghost(2) + +! Write and info message about what we are doing. + call CCTK_INFO ( 'Integrating equatorial circumference' ) + +! Setup the integration array for the equatorial circumference +! calculation. Note this is still a 2d-array, so we have to +! do the sum only for the correct 1d-slice of the array. + int_tmp = phi_sym_factor * phiweights * dlphi + +! If we are on a processor that contains part of the correct 1-d slice +! sum up the appropriate elemets, otherwise this processor will +! contribute zero. + if ( ( ithl+lbnd(1) .le. githeta ) .and. ( githeta .le. ithr+lbnd(1) ) ) then + itha = githeta - lbnd(1) + eq_circ_loc = sum ( int_tmp(itha,jphl:jphr) ) + else + eq_circ_loc = zero + end if + +! Then reduce the results over all processors and send the result to +! all processors. + call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & + eq_circ_loc, eq_circ, CCTK_VARIABLE_REAL ) + +! Get the upper and lower bounds for the circumference arrays. + call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd1, 'ehfinder::eh_circ_eq2' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get lower bounds for area array' ) + end if + call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd1, 'ehfinder::eh_circ_eq2' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get upper bounds for area array' ) + end if + call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh1, 'ehfinder::eh_circ_eq2' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for area array' ) + end if + +! If we are on the correct processor store the result. + if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then + k = sn - lbnd1(1) + eh_circ_eq2(k,l) = eq_circ + end if + +! Write an info message with the calculated circunference. + write(info_message,'(a35,f10.5)') 'Horizon equatorial circumference = ', & + eq_circ + call CCTK_INFO(info_message) + +! Now do the same for the polar circomference. + call CCTK_INFO ( 'Integrating polar circumference' ) + +! Setup the integration array for the polar circumference +! calculation. Note this is still a 2d-array, so we have to +! do the sum only for the correct 1d-slice of the array. + int_tmp = theta_sym_factor * thetaweights * dltheta + +! If we are on a processor that contains part of the correct 1-d slice +! sum up the appropriate elemets, otherwise this processor will +! contribute zero. + if ( ( jphl+lbnd(2) .le. gjphi ) .and. ( gjphi .le. jphr+lbnd(2) ) ) then + jpha = gjphi - lbnd(2) + pol_circ_loc = sum ( int_tmp(ithl:ithr,jpha) ) + else + pol_circ_loc = zero + end if + +! Then reduce the results over all processors and send the result to +! all processors. + call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & + pol_circ_loc, pol_circ, CCTK_VARIABLE_REAL ) + + +! If we are on the correct processor store the result. + if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then + k = sn - lbnd1(1) + eh_circ_pol2(k,l) = pol_circ + end if + +! Write an info message with the calculated circumference. + write(info_message,'(a30,f10.5)') 'Horizon polar circumference = ',pol_circ + call CCTK_INFO(info_message) + +end subroutine EHFinder_IntegrateCircumference + + +! This routine sets up the weights for the Simpsons rule integration +! over the surface. +subroutine EHFinder_InitWeights(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, im, jm + + CCTK_INT, dimension(4) :: bbox + CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost + +! Find out the lower bounds of the distributed integration grid arrays. + call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get lower bounds for surface arrays' ) + end if + +! Find out the local size of the distributed integration grid arrays + call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) + end if + + weights = one + +! Initialise the weight grid array for the 2D Simpsons rule integration. +! To do this I need to figure out the global location of the given point. +! There are 3 cases in the one dimensional case. If the point is on the +! boundary the weight is 1/3. If it is at an even position the weight +! is 4/3 and if it is at an odd position the weight is 2/3. + + do j = 1, lsh(2) + +! This is first done in the phi direction. Meaning that all points with +! the same theta coordinate are set to the same weight. + if ( ( lbnd(2)+j .eq. 1 ) .or. ( lbnd(2)+j .eq. nphi ) ) then + weights(:,j) = onethird + phiweights(:,j) = onethird + else if ( mod(lbnd(2)+j,2) .eq. 0 ) then + weights(:,j) = fourthirds + phiweights(:,j) = fourthirds + else + weights(:,j) = twothirds + phiweights(:,j) = twothirds + end if + +! Then it is done in the theta direction with the one-dimensional +! weights beeing multiplied. + do i = 1, lsh(1) + if ( ( lbnd(1)+i .eq. 1 ) .or. ( lbnd(1)+i .eq. ntheta ) ) then + weights(i,j) = onethird * weights(i,j) + thetaweights(i,j) = onethird + else if ( mod(lbnd(1)+i,2) .eq. 0 ) then + weights(i,j) = fourthirds * weights(i,j) + thetaweights(i,j) = fourthirds + else + weights(i,j) = twothirds * weights(i,j) + thetaweights(i,j) = twothirds + end if + end do + end do + +! The end result is a 2D array with the weights in the following pattern. + +! 1/9 4/9 2/9 4/9 2/9 4/9 1/9 +! 4/9 16/9 8/9 16/9 8/9 16/9 4/9 +! 2/9 8/9 4/9 8/9 4/9 8/9 2/9 +! 4/9 16/9 8/9 16/9 8/9 16/9 4/9 +! 2/9 8/9 4/9 8/9 4/9 8/9 2/9 +! 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 + + +! This is the routine where the stored areas (in eh_area2) are copied +! into the trigger variable eh_area. +subroutine EHFinder_CopyArea(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + +! If finding of surface failed set area to zero. + 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 + + +! This is the routine where the stored centroids (in eh_centroid2_[xyz]) +! are copied into the trigger variable eh_centroid_[xyz]. +subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + +! 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 + + +! This is the routine where the stored circumferences (in eh_circ_[eq,pol]2) +! are copied into the trigger variable eh_circ_[eq,pol]. +subroutine EHFinder_CopyCircumference(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + 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 |