From dfdb43d790c3aa6f77dc14a301396ea89ed1be70 Mon Sep 17 00:00:00 2001 From: knarf Date: Tue, 12 Jan 2010 21:13:14 +0000 Subject: add forgotten file git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Hydro_Analysis/trunk@37 b6729ddc-ac74-4bd1-908c-9dc7244c52a1 --- src/Whisky_Separation.F90 | 164 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 164 insertions(+) create mode 100644 src/Whisky_Separation.F90 diff --git a/src/Whisky_Separation.F90 b/src/Whisky_Separation.F90 new file mode 100644 index 0000000..6dc434c --- /dev/null +++ b/src/Whisky_Separation.F90 @@ -0,0 +1,164 @@ +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + + /*@@ + @routine Hydro_Analysis_FindSeparation + @date Thu May 20 12:35:20 2004 + @author Ian Hawke + @desc + Finds the separation (in coordinate and proper distance) + between the NS. Equal mass symmetry is assumed so it is + actually the distance between the origin and the location + of maximum density. This is along a straight line, not a + geodesic, so still not perfect. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine Hydro_Analysis_FindSeparation(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL :: separation_dx, separation_dy, separation_dz + CCTK_REAL, dimension(:), allocatable :: separation_x, separation_y, & + separation_z + CCTK_REAL, dimension(:), allocatable :: s_gxx, s_gxy, s_gxz, & + s_gyy, s_gyz, s_gzz + CCTK_INT :: ierr, i + CCTK_INT :: param_table_handle, interp_handle, coord_system_handle + CCTK_INT :: vindex + + CCTK_POINTER, dimension(3) :: interp_coords + CCTK_INT, dimension(6) :: in_array_indices + CCTK_POINTER, dimension(6) :: out_arrays + CCTK_INT, dimension(6) :: out_array_type_codes + +!!$ Proper separation requires interpolation + + allocate(separation_x(Hydro_Analysis_rho_max_origin_distance_npoints), & + separation_y(Hydro_Analysis_rho_max_origin_distance_npoints), & + separation_z(Hydro_Analysis_rho_max_origin_distance_npoints), STAT=ierr) + + if (ierr .ne. 0) then + call CCTK_WARN(0, "Failed to allocate separation coordinate arrays") + end if + + separation_dx = Hydro_Analysis_rho_max_loc(1) / dble(Hydro_Analysis_rho_max_origin_distance_npoints) + separation_dy = Hydro_Analysis_rho_max_loc(2) / dble(Hydro_Analysis_rho_max_origin_distance_npoints) + separation_dz = Hydro_Analysis_rho_max_loc(3) / dble(Hydro_Analysis_rho_max_origin_distance_npoints) + + do i = 1, Hydro_Analysis_rho_max_origin_distance_npoints + + separation_x(i) = i * separation_dx + separation_y(i) = i * separation_dy + separation_z(i) = i * separation_dz + + end do + + allocate(s_gxx(Hydro_Analysis_rho_max_origin_distance_npoints), & + s_gxy(Hydro_Analysis_rho_max_origin_distance_npoints), & + s_gxz(Hydro_Analysis_rho_max_origin_distance_npoints), & + s_gyy(Hydro_Analysis_rho_max_origin_distance_npoints), & + s_gyz(Hydro_Analysis_rho_max_origin_distance_npoints), & + s_gzz(Hydro_Analysis_rho_max_origin_distance_npoints), & + STAT=ierr) + + if (ierr .ne. 0) then + call CCTK_WARN(0, "Failed to allocate separation metric arrays") + end if + + param_table_handle = -1 + interp_handle = -1 + coord_system_handle = -1 + + call Util_TableCreateFromString (param_table_handle, "order = 2") + if (param_table_handle .lt. 0) then + call CCTK_WARN(0,"Cannot create parameter table for interpolator") + endif + + call CCTK_InterpHandle (interp_handle, "uniform cartesian") + if (interp_handle.lt.0) then + call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators (e.g. LocalInterp)?") + endif + + 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 (e.g. LocalInterp)?") + endif + +! fill in the input/output arrays for the interpolator + interp_coords(1) = CCTK_PointerTo(separation_x) + interp_coords(2) = CCTK_PointerTo(separation_y) + interp_coords(3) = CCTK_PointerTo(separation_z) + + call CCTK_VarIndex (vindex, "admbase::gxx") + in_array_indices(1) = vindex + call CCTK_VarIndex (vindex, "admbase::gyy") + in_array_indices(2) = vindex + call CCTK_VarIndex (vindex, "admbase::gzz") + in_array_indices(3) = vindex + call CCTK_VarIndex (vindex, "admbase::gxy") + in_array_indices(4) = vindex + call CCTK_VarIndex (vindex, "admbase::gxz") + in_array_indices(5) = vindex + call CCTK_VarIndex (vindex, "admbase::gyz") + in_array_indices(6) = vindex + + out_arrays(1) = CCTK_PointerTo(s_gxx) + out_arrays(2) = CCTK_PointerTo(s_gyy) + out_arrays(3) = CCTK_PointerTo(s_gzz) + out_arrays(4) = CCTK_PointerTo(s_gxy) + out_arrays(5) = CCTK_PointerTo(s_gxz) + out_arrays(6) = CCTK_PointerTo(s_gyz) + + out_array_type_codes = CCTK_VARIABLE_REAL + +! Interpolation. + + call CCTK_InterpGridArrays (ierr, cctkGH, 3, interp_handle,& + param_table_handle, coord_system_handle,& + Hydro_Analysis_rho_max_origin_distance_npoints, CCTK_VARIABLE_REAL, interp_coords,& + 6, in_array_indices,& + 6, out_array_type_codes, out_arrays) + if (ierr < 0) then + call CCTK_WARN (1, "interpolator call returned an error code"); + endif + +! release parameter table + call Util_TableDestroy (ierr, param_table_handle) + +! Integrate using trapezoidal rule. + + Hydro_Analysis_rho_max_origin_distance = 0.d0 + + do i=2, Hydro_Analysis_rho_max_origin_distance_npoints + + Hydro_Analysis_rho_max_origin_distance = Hydro_Analysis_rho_max_origin_distance + 0.5d0 & + *(sqrt(s_gxx(i-1)*separation_dx**2 + & + s_gyy(i-1)*separation_dy**2 + & + s_gzz(i-1)*separation_dz**2 & + + 2.d0*(s_gxy(i-1)*separation_dx*separation_dy + & + s_gxz(i-1)*separation_dx*separation_dz + & + s_gyz(i-1)*separation_dy*separation_dz)) & + + sqrt(s_gxx(i )*separation_dx**2 + & + s_gyy(i )*separation_dy**2 + & + s_gzz(i )*separation_dz**2 & + + 2.d0*(s_gxy(i )*separation_dx*separation_dy + & + s_gxz(i )*separation_dx*separation_dz + & + s_gyz(i )*separation_dy*separation_dz))) + + end do + +end subroutine Hydro_Analysis_FindSeparation + -- cgit v1.2.3