From 43d765ed9db6ebd283836cd9b58ce01daa05147f Mon Sep 17 00:00:00 2001 From: diener Date: Mon, 26 Jan 2009 23:54:25 +0000 Subject: Use modules to compile on the SiCortex machine (from Steve Brandt). git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@202 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_IsoSurface_optimized.F90 | 303 ++++++++++++++++++---------------- 1 file changed, 158 insertions(+), 145 deletions(-) diff --git a/src/EHFinder_IsoSurface_optimized.F90 b/src/EHFinder_IsoSurface_optimized.F90 index 90ba5d6..18641a9 100755 --- a/src/EHFinder_IsoSurface_optimized.F90 +++ b/src/EHFinder_IsoSurface_optimized.F90 @@ -6,150 +6,8 @@ #include "cctk_Arguments.h" #include "cctk_Functions.h" -subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS) - - use EHFinder_mod - use IsoSurface_mod - - implicit none - - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - CCTK_REAL, dimension(:), pointer :: x_coord, y_coord, z_coord - CCTK_INT, dimension(:), pointer :: t_index - CCTK_INT :: l,contor,p,nsize - -! Find index ranges for points excluding ghost zones and symmetry points for -! the 3D grid functions. -#include "include/physical_part.h" - - call CCTK_INFO ('Entered IsoSurface') - - nx = cctk_lsh(1); ny = cctk_lsh(2); nz = cctk_lsh(3) - - nsize=1000 - allocate(x_coord(nsize),y_coord(nsize),z_coord(nsize),t_index(5*nsize)) - - do p = 1, eh_number_level_sets - call triangularization(l,contor,nx,ny,nz,x_coord,y_coord,z_coord, & - t_index,f(:,:,:,p),x(:,1,1),y(1,:,1),z(1,1,:)) - - print*, 'contor = ', t_index(1:contor) - print*,'The surface consists of', contor/3,'triangles and',l,'coordinates' - - call interpolation_area(cctkGH, l, x_coord,y_coord,z_coord, contor, & - t_index, eh_area(1,p) ) - - print*,'The area is ', eh_area(1,p) - end do - - - -end subroutine EHFinder_IsoSurface - - -subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,area) - - use EHFinder_mod - implicit none - DECLARE_CCTK_PARAMETERS -! DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - CCTK_POINTER, intent(IN) :: cctkGH - CCTK_INT, intent(IN) :: npoints, l - CCTK_REAL, dimension(:), pointer :: x_coord,y_coord,z_coord - CCTK_INT, dimension(:), pointer :: t_index - CCTK_REAL, intent(OUT) :: area - character(len=200) :: area_interp - character(len=128) :: warn_message - CCTK_INT :: area_interp_len, local_contor - character(len=7) :: area_order - CCTK_INT :: interp_handle, table_handle, coord_system_handle - CCTK_POINTER, dimension(3) :: interp_coords - CCTK_POINTER, dimension(6) :: out_array - CCTK_INT, dimension(6) :: in_array - CCTK_REAL, dimension(:), allocatable :: metric_xx, metric_yy, metric_zz, metric_xy, metric_yz, metric_xz - - CCTK_INT, dimension(6), parameter :: out_types = (/ CCTK_VARIABLE_REAL, & - CCTK_VARIABLE_REAL, & - CCTK_VARIABLE_REAL, & - CCTK_VARIABLE_REAL, & - CCTK_VARIABLE_REAL, & - CCTK_VARIABLE_REAL /) - ! Convert the interpolator parameter into a fortran string. - call CCTK_FortranString ( area_interp_len, area_interpolator, & - area_interp ) - -!allocation of metric arrays - allocate(metric_xx(npoints), metric_yy(npoints), metric_zz(npoints), metric_xy(npoints), metric_xz(npoints), metric_yz(npoints)) - -! 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 = 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 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 = trim(warn_message)//'Forgot to activate an implementation ' - warn_message = trim(warn_message)//'providing coordinates?' - call CCTK_WARN( 0, trim(warn_message) ) - endif - -! Get the pointers to the interpolation points. - interp_coords(1) = CCTK_PointerTo(x_coord) - interp_coords(2) = CCTK_PointerTo(y_coord) - interp_coords(3) = CCTK_PointerTo(z_coord) - -! Get the pointers to the interpolation return arrays. - out_array(1) = CCTK_PointerTo(metric_xx) - out_array(2) = CCTK_PointerTo(metric_xy) - out_array(3) = CCTK_PointerTo(metric_xz) - out_array(4) = CCTK_PointerTo(metric_yy) - out_array(5) = CCTK_PointerTo(metric_yz) - out_array(6) = CCTK_PointerTo(metric_zz) -! out_array(7) = CCTK_PointerTo(psii) -! 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' ) - - call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, & - table_handle, coord_system_handle, & - npoints, CCTK_VARIABLE_REAL, & - interp_coords, 6, in_array(1:6), & - 6, out_types(1:6), out_array(1:6) ) - - call area_computation(x_coord, y_coord, z_coord, t_index, & - metric_xx, metric_xy, metric_xz, & - metric_yy, metric_yz, metric_zz, & - npoints, l, area ) - - deallocate(metric_xx, metric_yy, metric_zz, metric_xz, metric_yz, metric_xy) - - end subroutine interpolation_area - + module areac + contains subroutine area_computation(x_coord, y_coord, z_coord, t_index, & g11, g12, g13, & g22, g23, g33, & @@ -203,8 +61,11 @@ subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,a area=area+ds end do end subroutine area_computation + end module - + + module triangular + contains subroutine triangularization(n_coords,n_triangles,Nx,Ny,Nz,x_coord,y_coord,z_coord,t_index,grid_value,grid_x,grid_y,grid_z) !subroutine triangularization(l,contor,nx,ny,nz,x_coord,y_coord,z_coord,t_index,f(:,:,:,1),x(:,1,1),y(1,:,1),z(1,1,:)) @@ -822,3 +683,155 @@ integer, dimension (16,16):: b9(16,16)=reshape& close(12) print*,'t_index = ',t_index(1:n_triangles) end subroutine triangularization + end module + +module interpa +contains +subroutine interpolation_area(cctkGH,npoints,x_coord,y_coord,z_coord,l,t_index,area) + + use areac + use EHFinder_mod + implicit none + DECLARE_CCTK_PARAMETERS +! DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_POINTER, intent(IN) :: cctkGH + CCTK_INT, intent(IN) :: npoints, l + CCTK_REAL, dimension(:), pointer :: x_coord,y_coord,z_coord + CCTK_INT, dimension(:), pointer :: t_index + CCTK_REAL, intent(OUT) :: area + character(len=200) :: area_interp + character(len=128) :: warn_message + CCTK_INT :: area_interp_len, local_contor + character(len=7) :: area_order + CCTK_INT :: interp_handle, table_handle, coord_system_handle + CCTK_POINTER, dimension(3) :: interp_coords + CCTK_POINTER, dimension(6) :: out_array + CCTK_INT, dimension(6) :: in_array + CCTK_REAL, dimension(:), allocatable :: metric_xx, metric_yy, metric_zz, metric_xy, metric_yz, metric_xz + + CCTK_INT, dimension(6), parameter :: out_types = (/ CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL /) + ! Convert the interpolator parameter into a fortran string. + call CCTK_FortranString ( area_interp_len, area_interpolator, & + area_interp ) + +!allocation of metric arrays + allocate(metric_xx(npoints), metric_yy(npoints), metric_zz(npoints), metric_xy(npoints), metric_xz(npoints), metric_yz(npoints)) + +! 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 = 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 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 = trim(warn_message)//'Forgot to activate an implementation ' + warn_message = trim(warn_message)//'providing coordinates?' + call CCTK_WARN( 0, trim(warn_message) ) + endif + +! Get the pointers to the interpolation points. + interp_coords(1) = CCTK_PointerTo(x_coord) + interp_coords(2) = CCTK_PointerTo(y_coord) + interp_coords(3) = CCTK_PointerTo(z_coord) + +! Get the pointers to the interpolation return arrays. + out_array(1) = CCTK_PointerTo(metric_xx) + out_array(2) = CCTK_PointerTo(metric_xy) + out_array(3) = CCTK_PointerTo(metric_xz) + out_array(4) = CCTK_PointerTo(metric_yy) + out_array(5) = CCTK_PointerTo(metric_yz) + out_array(6) = CCTK_PointerTo(metric_zz) +! out_array(7) = CCTK_PointerTo(psii) +! 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' ) + + call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + npoints, CCTK_VARIABLE_REAL, & + interp_coords, 6, in_array(1:6), & + 6, out_types(1:6), out_array(1:6) ) + + call area_computation(x_coord, y_coord, z_coord, t_index, & + metric_xx, metric_xy, metric_xz, & + metric_yy, metric_yz, metric_zz, & + npoints, l, area ) + + deallocate(metric_xx, metric_yy, metric_zz, metric_xz, metric_yz, metric_xy) + + end subroutine interpolation_area + end module + +subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS) + + use interpa + use EHFinder_mod + use IsoSurface_mod + use triangular + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL, dimension(:), pointer :: x_coord, y_coord, z_coord + CCTK_INT, dimension(:), pointer :: t_index + CCTK_INT :: l,contor,p,nsize + +! Find index ranges for points excluding ghost zones and symmetry points for +! the 3D grid functions. +#include "include/physical_part.h" + + call CCTK_INFO ('Entered IsoSurface') + + nx = cctk_lsh(1); ny = cctk_lsh(2); nz = cctk_lsh(3) + + nsize=1000 + allocate(x_coord(nsize),y_coord(nsize),z_coord(nsize),t_index(5*nsize)) + + do p = 1, eh_number_level_sets + call triangularization(l,contor,nx,ny,nz,x_coord,y_coord,z_coord, & + t_index,f(:,:,:,p),x(:,1,1),y(1,:,1),z(1,1,:)) + + print*, 'contor = ', t_index(1:contor) + print*,'The surface consists of', contor/3,'triangles and',l,'coordinates' + + call interpolation_area(cctkGH, l, x_coord,y_coord,z_coord, contor, & + t_index, eh_area(1,p) ) + + print*,'The area is ', eh_area(1,p) + end do + + + +end subroutine EHFinder_IsoSurface + + -- cgit v1.2.3