diff --git a/doc/USER_MANUAL/bibliography.bib b/doc/USER_MANUAL/bibliography.bib index 5ad8cb198..7a98def6a 100644 --- a/doc/USER_MANUAL/bibliography.bib +++ b/doc/USER_MANUAL/bibliography.bib @@ -13621,6 +13621,14 @@ @Article{MeViSa06 Volume = {166} } +@article{Meschede2011, + author = {Meschede, M.A. and C.L. Myhrvold and J. Tromp}, + journal = gji, + pages = {529-537}, + title = {Antipodal focusing of seismic waves due to large meteorite impacts on Earth}, + volume = {187}, + year = {2011}} + @Article{Merkel07, Title = {Deformation of ({M}g,{F}e){S}i{O}3 post-perovskite and {D}'' anisotropy}, Author = {Merkel, S. and A. K. McNamara and A. Kubo and S. Speziale and L. Miyagi and Y. Meng and T. S. Duffy and H. R. Wenk}, diff --git a/doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf b/doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf index 2e512a738..53b6dc72a 100644 Binary files a/doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf and b/doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf differ diff --git a/setup/constants.h.in b/setup/constants.h.in index ab1e216cf..09b34e89a 100644 --- a/setup/constants.h.in +++ b/setup/constants.h.in @@ -220,6 +220,9 @@ ! (brute-force search will be slower but even more accurate) logical, parameter :: DO_BRUTE_FORCE_POINT_SEARCH = .false. +! minimum number of stations for using hash table in duplets search + integer, parameter :: DUPLETS_NREC_MINIMUM_FOR_HASH = 10000 + !!----------------------------------------------------------- !! !! CPML perfectly matched absorbing layers diff --git a/src/shared/read_topo_bathy_file.f90 b/src/shared/read_topo_bathy_file.f90 index 9c8e84bac..ffe5c373b 100644 --- a/src/shared/read_topo_bathy_file.f90 +++ b/src/shared/read_topo_bathy_file.f90 @@ -135,7 +135,11 @@ subroutine get_topo_elevation_free(x_target,y_target,target_elevation,target_dis ! get approximate topography elevation at source long/lat coordinates - use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NGLLSQUARE,HUGEVAL,TINYVAL,MIDX,MIDY,MIDZ,USE_DISTANCE_CRITERION_TOPO + use constants, only: myrank,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NGLLSQUARE, & + HUGEVAL,TINYVAL,MIDX,MIDY,MIDZ, & + USE_DISTANCE_CRITERION_TOPO + + use shared_parameters, only: free_surface_xy_midpoints,free_surface_typical_size,free_surface_has_typical_size implicit none @@ -162,12 +166,11 @@ subroutine get_topo_elevation_free(x_target,y_target,target_elevation,target_dis real(kind=CUSTOM_REAL) :: distmin,dist real(kind=CUSTOM_REAL) :: dist_node_min,norm - integer :: iface,i,j,k,ispec,iglob,igll,jgll,kgll,ijk + integer :: iface,i,j,k,ispec,iglob,igll,jgll,kgll,ijk,ier integer :: iselected,jselected,iface_selected integer :: inode,iadjust,jadjust integer :: ilocmin(1) - real(kind=CUSTOM_REAL) :: typical_size logical :: located_target ! initialize @@ -178,11 +181,30 @@ subroutine get_topo_elevation_free(x_target,y_target,target_elevation,target_dis ! computes typical size of elements at the surface (uses first element for estimation) if (USE_DISTANCE_CRITERION_TOPO) then - ispec = free_surface_ispec(1) - typical_size = (xstore(ibool(1,1,1,ispec)) - xstore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 & - + (ystore(ibool(1,1,1,ispec)) - ystore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 - ! use 10 times the distance as a criterion for point detection - typical_size = 10.0_CUSTOM_REAL * typical_size + ! do only once + if (.not. free_surface_has_typical_size) then + ispec = free_surface_ispec(1) + free_surface_typical_size = (xstore(ibool(1,1,1,ispec)) - xstore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 & + + (ystore(ibool(1,1,1,ispec)) - ystore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 + ! use 10 times the distance as a criterion for point detection + free_surface_typical_size = 10.0_CUSTOM_REAL * free_surface_typical_size + + ! prepares midpoints coordinates + allocate(free_surface_xy_midpoints(2,num_free_surface_faces),stat=ier) + if (ier /= 0 ) call exit_MPI(myrank,'Error allocating array free_surface_xy_midpoints') + free_surface_xy_midpoints(:,:) = 0.0_CUSTOM_REAL + + ! store x/y coordinates of center point + do iface = 1,num_free_surface_faces + ispec = free_surface_ispec(iface) + iglob = ibool(MIDX,MIDY,MIDZ,ispec) + free_surface_xy_midpoints(1,iface) = xstore(iglob) + free_surface_xy_midpoints(2,iface) = ystore(iglob) + enddo + + ! done setup + free_surface_has_typical_size = .true. + endif endif ! flag to check that we located at least one target element @@ -196,15 +218,17 @@ subroutine get_topo_elevation_free(x_target,y_target,target_elevation,target_dis ! loops over all free surface faces do iface = 1,num_free_surface_faces - ispec = free_surface_ispec(iface) ! exclude elements that are too far from target if (USE_DISTANCE_CRITERION_TOPO) then - iglob = ibool(MIDX,MIDY,MIDZ,ispec) - dist = (x_target - xstore(iglob))**2 + (y_target - ystore(iglob))**2 - if (dist > typical_size) cycle + dist = (x_target - free_surface_xy_midpoints(1,iface))*(x_target - free_surface_xy_midpoints(1,iface)) & + + (y_target - free_surface_xy_midpoints(2,iface))*(y_target - free_surface_xy_midpoints(2,iface)) + + if (dist > free_surface_typical_size) cycle endif + ispec = free_surface_ispec(iface) + ! loop only on points inside the element ! exclude edges to ensure this point is not shared with other elements do j = 2,NGLLY - 1 @@ -217,8 +241,8 @@ subroutine get_topo_elevation_free(x_target,y_target,target_elevation,target_dis iglob = ibool(igll,jgll,kgll,ispec) ! distance (squared) to target - dist = ( x_target - xstore(iglob) )**2 + & - ( y_target - ystore(iglob) )**2 + dist = (x_target - xstore(iglob))*(x_target - xstore(iglob)) & + + (y_target - ystore(iglob))*(y_target - ystore(iglob)) ! keep this point if it is closer to the receiver if (dist < distmin) then @@ -284,17 +308,22 @@ subroutine get_topo_elevation_free(x_target,y_target,target_elevation,target_dis ! stores node infos inode = inode + 1 elevation_node(inode) = zstore(iglob) - dist_node(inode) = sqrt( (x_target - xstore(iglob))**2 + (y_target - ystore(iglob))**2) + + ! distance squared + dist_node(inode) = (x_target - xstore(iglob))*(x_target - xstore(iglob)) & + + (y_target - ystore(iglob))*(y_target - ystore(iglob)) enddo enddo ! weighted elevation - dist = sum( dist_node(:) ) + dist = dist_node(1) + dist_node(2) + dist_node(3) + dist_node(4) if (dist < distmin) then ! sets new minimum distance (of all 4 closest nodes) distmin = dist - target_distmin = distmin + + ! outputs sum of squared distances (used to evaluate closest mesh locations) + target_distmin = dist ! interpolates elevation if (dist > TINYVAL) then @@ -309,7 +338,7 @@ subroutine get_topo_elevation_free(x_target,y_target,target_elevation,target_dis ! gets minimum distance value & index ilocmin = minloc(dist_node(:)) - dist_node_min = dist_node(ilocmin(1)) + dist_node_min = sqrt(dist_node(ilocmin(1))) ! checks if a node is almost exactly at target location if (dist_node_min < TINYVAL) then @@ -320,19 +349,20 @@ subroutine get_topo_elevation_free(x_target,y_target,target_elevation,target_dis ! (closer points have higher weight) do k = 1,4 if (dist_node(k) > 0.0_CUSTOM_REAL) then - weight(k) = 1.0_CUSTOM_REAL / dist_node(k) + ! needs exact distance for interpolation weight + weight(k) = 1.0_CUSTOM_REAL / sqrt(dist_node(k)) else weight(k) = 1.e10 ! very large weight for point on target, will dominate interpolation endif enddo ! normalize weights: w_i = w_i / sum(w_i) - norm = sum(weight(:)) + norm = weight(1) + weight(2) + weight(3) + weight(4) if (norm > TINYVAL) then weight(:) = weight(:) / norm else ! all 4 weights almost zero, meaning distances are all very large; uses equal weighting - weight(:) = 1.0 / 4.0 + weight(:) = 0.25_CUSTOM_REAL ! 1.0 / 4.0 endif ! interpolation target_elevation = weight(1)*elevation_node(1) & @@ -362,7 +392,11 @@ subroutine get_topo_elevation_free_closest(x_target,y_target,target_elevation,ta ! get approximate topography elevation at long/lat coordinates from closest point - use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NGLLSQUARE,HUGEVAL,MIDX,MIDY,MIDZ,USE_DISTANCE_CRITERION_TOPO + use constants, only: myrank,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NGLLSQUARE, & + HUGEVAL,MIDX,MIDY,MIDZ, & + USE_DISTANCE_CRITERION_TOPO + + use shared_parameters, only: free_surface_xy_midpoints,free_surface_typical_size,free_surface_has_typical_size implicit none @@ -386,9 +420,7 @@ subroutine get_topo_elevation_free_closest(x_target,y_target,target_elevation,ta ! local parameters real(kind=CUSTOM_REAL) :: distmin,dist - integer :: iface,i,ispec,iglob,igll,jgll,kgll - - real(kind=CUSTOM_REAL) :: typical_size + integer :: iface,i,ispec,iglob,igll,jgll,kgll,ier logical :: located_target ! initialize @@ -399,11 +431,30 @@ subroutine get_topo_elevation_free_closest(x_target,y_target,target_elevation,ta ! computes typical size of elements at the surface (uses first element for estimation) if (USE_DISTANCE_CRITERION_TOPO) then - ispec = free_surface_ispec(1) - typical_size = (xstore(ibool(1,1,1,ispec)) - xstore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 & - + (ystore(ibool(1,1,1,ispec)) - ystore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 - ! use 10 times the distance as a criterion for point detection - typical_size = 10.0_CUSTOM_REAL * typical_size + ! do only once + if (.not. free_surface_has_typical_size) then + ispec = free_surface_ispec(1) + free_surface_typical_size = (xstore(ibool(1,1,1,ispec)) - xstore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 & + + (ystore(ibool(1,1,1,ispec)) - ystore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 + ! use 10 times the distance as a criterion for point detection + free_surface_typical_size = 10.0_CUSTOM_REAL * free_surface_typical_size + + ! prepares midpoints coordinates + allocate(free_surface_xy_midpoints(2,num_free_surface_faces),stat=ier) + if (ier /= 0 ) call exit_MPI(myrank,'Error allocating array free_surface_xy_midpoints') + free_surface_xy_midpoints(:,:) = 0.0_CUSTOM_REAL + + ! store x/y coordinates of center point + do iface = 1,num_free_surface_faces + ispec = free_surface_ispec(iface) + iglob = ibool(MIDX,MIDY,MIDZ,ispec) + free_surface_xy_midpoints(1,iface) = xstore(iglob) + free_surface_xy_midpoints(2,iface) = ystore(iglob) + enddo + + ! done setup + free_surface_has_typical_size = .true. + endif endif ! flag to check that we located at least one target element @@ -414,15 +465,17 @@ subroutine get_topo_elevation_free_closest(x_target,y_target,target_elevation,ta ! loops over all free surface faces do iface = 1,num_free_surface_faces - ispec = free_surface_ispec(iface) ! excludes elements that are too far from target if (USE_DISTANCE_CRITERION_TOPO) then - iglob = ibool(MIDX,MIDY,MIDZ,ispec) - dist = (x_target - xstore(iglob))**2 + (y_target - ystore(iglob))**2 - if (dist > typical_size) cycle + dist = (x_target - free_surface_xy_midpoints(1,iface))*(x_target - free_surface_xy_midpoints(1,iface)) & + + (y_target - free_surface_xy_midpoints(2,iface))*(y_target - free_surface_xy_midpoints(2,iface)) + + if (dist > free_surface_typical_size) cycle endif + ispec = free_surface_ispec(iface) + ! loop only on points inside the element do i = 1,NGLLSQUARE igll = free_surface_ijk(1,i,iface) @@ -432,8 +485,8 @@ subroutine get_topo_elevation_free_closest(x_target,y_target,target_elevation,ta iglob = ibool(igll,jgll,kgll,ispec) ! distance (squared) to target - dist = ( x_target - xstore(iglob) )**2 + & - ( y_target - ystore(iglob) )**2 + dist = (x_target - xstore(iglob))*(x_target - xstore(iglob)) & + + (y_target - ystore(iglob))*(y_target - ystore(iglob)) ! keep this point if it is closer to the receiver if (dist < distmin) then @@ -455,7 +508,8 @@ subroutine get_topo_elevation_free_closest(x_target,y_target,target_elevation,ta iglob = ibool(1,1,1,ispec) ! elevation (given in z - coordinate) target_elevation = zstore(iglob) - target_distmin = ( x_target - xstore(iglob) )**2 + ( y_target - ystore(iglob) )**2 + target_distmin = (x_target - xstore(iglob))*(x_target - xstore(iglob)) & + + (y_target - ystore(iglob))*(y_target - ystore(iglob)) located_target = .true. endif diff --git a/src/shared/shared_par.F90 b/src/shared/shared_par.F90 index 532b01bd8..36ab72cde 100644 --- a/src/shared/shared_par.F90 +++ b/src/shared/shared_par.F90 @@ -213,6 +213,7 @@ end module shared_input_parameters module shared_compute_parameters ! parameters to be computed based upon parameters above read from file + use constants, only: CUSTOM_REAL implicit none @@ -236,6 +237,13 @@ module shared_compute_parameters ! fault rupture simulation logical :: FAULT_SIMULATION = .false. + ! free surface + ! for elevation search: x/y coordinates of free surface element midpoints + real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: free_surface_xy_midpoints + real(kind=CUSTOM_REAL) :: free_surface_typical_size + ! flag to calculate typical size only once + logical :: free_surface_has_typical_size = .false. + end module shared_compute_parameters ! diff --git a/src/specfem3D/get_elevation.f90 b/src/specfem3D/get_elevation.f90 index dd5772d06..cbbb1ed08 100644 --- a/src/specfem3D/get_elevation.f90 +++ b/src/specfem3D/get_elevation.f90 @@ -120,6 +120,7 @@ subroutine get_elevation_and_z_coordinate_all(npoints,plon,plat,pbur,utm_x,utm_y elevation(ipoin) = loc_ele elevation_distmin(ipoin) = loc_distmin enddo + ! MPI communications to determine the best slice call gather_all_dp(elevation,npoints,elevation_all,npoints,NPROC) call gather_all_dp(elevation_distmin,npoints,elevation_distmin_all,npoints,NPROC) diff --git a/src/specfem3D/locate_point.f90 b/src/specfem3D/locate_point.f90 index a55a457f4..203ee0c81 100644 --- a/src/specfem3D/locate_point.f90 +++ b/src/specfem3D/locate_point.f90 @@ -38,7 +38,8 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & use specfem_par, only: ibool,myrank,NSPEC_AB,NGLOB_AB, & xstore,ystore,zstore, & ispec_is_surface_external_mesh,iglob_is_surface_external_mesh, & - xyz_midpoints + xyz_midpoints, & + ACOUSTIC_SIMULATION,ELASTIC_SIMULATION use specfem_par_acoustic, only: ispec_is_acoustic use specfem_par_elastic, only: ispec_is_elastic @@ -72,10 +73,13 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & integer :: imin,imax,jmin,jmax,kmin,kmax !! DK DK dec 2017: also loop on all the elements in contact with the initial guess element to improve accuracy of estimate - logical, dimension(:),allocatable :: flag_topological ! making array allocatable, otherwise will crash for large meshes + logical, dimension(:),allocatable, save :: flag_topological ! making array allocatable, otherwise will crash for large meshes integer :: number_of_mesh_elements_for_the_initial_guess logical :: use_adjacent_elements_search + ! flag to allocate topological array only once + logical, save :: has_flag_topological = .false. + ! dynamic array ! integer, dimension(:), allocatable :: array_of_all_elements_of_ispec_selected ! logical, parameter :: USE_SINGLE_PASS = .false. @@ -93,6 +97,14 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & logical :: use_brute_force_search logical :: is_better_location + logical :: do_include_slice + + ! flag to check mesh dimensions only once + logical, save :: has_slice_dimensions = .false. + double precision, save :: x_min_slice,x_max_slice + double precision, save :: y_min_slice,y_max_slice + double precision, save :: z_min_slice,z_max_slice + !debug !print *,'locate point in mesh: ',x_target, y_target, z_target @@ -114,6 +126,53 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & ! set distance to huge initial value distmin_squared = HUGEVAL + ! gets local slice dimension + if (.not. has_slice_dimensions) then + ! mesh slice dimension of current process + x_min_slice = minval(xstore) + x_max_slice = maxval(xstore) + y_min_slice = minval(ystore) + y_max_slice = maxval(ystore) + z_min_slice = minval(zstore) + z_max_slice = maxval(zstore) + ! only do this once + has_slice_dimensions = .true. + endif + + ! checks if target location is in this mesh slice + do_include_slice = .false. + + if (POINT_CAN_BE_BURIED .eqv. .false.) then + ! only checks if x/y position is in slice + if ((x_target < x_max_slice + elemsize_max_glob .and. x_target > x_min_slice - elemsize_max_glob) & + .and. (y_target < y_max_slice + elemsize_max_glob .and. y_target > y_min_slice - elemsize_max_glob)) then + do_include_slice = .true. + endif + else + ! checks x/y/z position + if ((x_target < x_max_slice + elemsize_max_glob .and. x_target > x_min_slice - elemsize_max_glob) & + .and. (y_target < y_max_slice + elemsize_max_glob .and. y_target > y_min_slice - elemsize_max_glob) & + .and. (z_target < z_max_slice + elemsize_max_glob .and. z_target > z_min_slice - elemsize_max_glob)) then + do_include_slice = .true. + endif + endif + + if (.not. do_include_slice) then + ! point outside of this mesh slice + ispec_selected = 0 + domain = 0 + + x_found = -HUGEVAL + y_found = -HUGEVAL + z_found = -HUGEVAL + + ! store final distance squared between asked and found + final_distance_squared = HUGEVAL + + ! done search in this slice + return + endif + ! point search type if (DO_BRUTE_FORCE_POINT_SEARCH .or. (POINT_CAN_BE_BURIED .eqv. .false.)) then use_brute_force_search = .true. @@ -179,7 +238,9 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & endif ! skip buried elements in case station must lie on surface - if (.not. POINT_CAN_BE_BURIED .and. .not. ispec_is_surface_external_mesh(ispec)) cycle + if (.not. POINT_CAN_BE_BURIED) then + if (.not. ispec_is_surface_external_mesh(ispec)) cycle + endif ! distance to element midpoint dist_squared = (x_target - xyz_midpoints(1,ispec))*(x_target - xyz_midpoints(1,ispec)) & @@ -193,10 +254,13 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & do k = kmin, kmax do j = jmin, jmax do i = imin, imax + iglob = ibool(i,j,k,ispec) ! skip inner GLL points in case station must lie on surface - if (.not. POINT_CAN_BE_BURIED .and. .not. iglob_is_surface_external_mesh(iglob)) cycle + if (.not. POINT_CAN_BE_BURIED) then + if (.not. iglob_is_surface_external_mesh(iglob)) cycle + endif x = dble(xstore(iglob)) y = dble(ystore(iglob)) @@ -277,6 +341,18 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & if (use_adjacent_elements_search) then !! DK DK dec 2017: also loop on all the elements in contact with the initial guess element to improve accuracy of estimate + ! allocates arrays + ! we only do this the first time running through this, and keep the arrays to speed this for the next point location + if (.not. has_flag_topological) then + ! topological flags to find neighboring elements + allocate(flag_topological(NGLOB_AB),stat=ier) + if (ier /= 0) stop 'Error allocating flag_topological array' + flag_topological(:) = .false. + + ! allocate only once + has_flag_topological = .true. + endif + if (use_brute_force_search) then ! brute-force search always loops over whole mesh slice num_elem_local = NSPEC_AB @@ -294,7 +370,7 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & kdtree_search_num_nodes = num_elem_local ! debug - !print *,' total number of search elements: ',num_elem_local,POINT_CAN_BE_BURIED,r_search,elemsize_max_glob + !print *,'debug: total number of search elements: ',num_elem_local,POINT_CAN_BE_BURIED,r_search,elemsize_max_glob ! allocates search array if (kdtree_search_num_nodes > 0) then @@ -320,13 +396,11 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & ! dummy search in this slice num_elem_local = 1 - kdtree_search_num_nodes = num_elem_local + kdtree_search_num_nodes = 1 endif endif ! flagging corners - allocate(flag_topological(NGLOB_AB),stat=ier) - if (ier /= 0) stop 'Error allocating flag_topological array' flag_topological(:) = .false. ! mark the eight corners of the initial guess element @@ -368,11 +442,13 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & ! this element is in contact with the initial guess number_of_mesh_elements_for_the_initial_guess = number_of_mesh_elements_for_the_initial_guess + 1 ! check - if (number_of_mesh_elements_for_the_initial_guess > 100) stop 'Error must increase array size in locate_point.f90' + if (number_of_mesh_elements_for_the_initial_guess > 100) & + stop 'Error must increase array size in locate_point.f90' array_of_all_elements_of_ispec_selected(number_of_mesh_elements_for_the_initial_guess) = ispec - ! let us not count it more than once, it may have a full edge in contact with it and would then be counted twice + ! let us not count it more than once, + ! it may have a full edge in contact with it and would then be counted twice goto 707 endif @@ -409,7 +485,8 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & ! this element is in contact with the initial guess number_of_mesh_elements_for_the_initial_guess = number_of_mesh_elements_for_the_initial_guess + 1 - ! let us not count it more than once, it may have a full edge in contact with it and would then be counted twice + ! let us not count it more than once, + ! it may have a full edge in contact with it and would then be counted twice goto 700 endif @@ -452,7 +529,8 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & number_of_mesh_elements_for_the_initial_guess = number_of_mesh_elements_for_the_initial_guess + 1 array_of_all_elements_of_ispec_selected(number_of_mesh_elements_for_the_initial_guess) = ispec - ! let us not count it more than once, it may have a full edge in contact with it and would then be counted twice + ! let us not count it more than once, + ! it may have a full edge in contact with it and would then be counted twice goto 800 endif @@ -476,7 +554,8 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & endif endif - deallocate(flag_topological) + ! we keep topological flags for next point locations to speed up routine, thus keep this here commented out. + ! deallocate(flag_topological) else ! no need for adjacent element, only loop within initial guess @@ -519,9 +598,11 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & is_better_location = .true. endif ! takes position if old position is in acoustic element and new one in elastic one - ! prefers having station in elastic elements over acoustic at interfaces - if (ispec_is_acoustic(ispec_selected) .and. ispec_is_elastic(ispec)) then - is_better_location = .true. + if (ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION) then + ! prefers having station in elastic elements over acoustic at interfaces + if (ispec_is_acoustic(ispec_selected) .and. ispec_is_elastic(ispec)) then + is_better_location = .true. + endif endif endif endif diff --git a/src/specfem3D/locate_receivers.f90 b/src/specfem3D/locate_receivers.f90 index c2ce13c14..cff3a40aa 100644 --- a/src/specfem3D/locate_receivers.f90 +++ b/src/specfem3D/locate_receivers.f90 @@ -92,6 +92,7 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) double precision, dimension(NDIM,NDIM,NREC_SUBSET_MAX) :: nu_subset integer, dimension(NREC_SUBSET_MAX) :: ispec_selected_rec_subset,idomain_subset integer :: nrec_subset_current_size,irec_in_this_subset,irec_already_done + integer :: num_output_info logical :: is_done_stations @@ -114,6 +115,10 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) write(IMAIN,*) ' (depth) becomes directly (z) coordinate' endif call flush_IMAIN() + + ! output frequency for large number of receivers + ! number to output about ~50 steps, rounds to the next multiple of 500 + num_output_info = max(500,int(ceiling(ceiling(nrec/50.0)/500.0)*500)) endif ! compute typical size of elements @@ -235,8 +240,9 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) ! user output progress if (myrank == 0 .and. nrec > 1000) then - if (mod(irec,500) == 0) then - write(IMAIN,*) ' located receivers ',irec,'out of',nrec + if (mod(irec,num_output_info) == 0) then + tCPU = wtime() - tstart + write(IMAIN,*) ' located receivers ',irec,'out of',nrec,' - elapsed time: ',sngl(tCPU),'s' call flush_IMAIN() endif endif @@ -398,17 +404,19 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) final_distance_max = maxval(final_distance(:)) ! display maximum error for all the receivers + write(IMAIN,*) write(IMAIN,*) 'maximum error in location of all the receivers: ',sngl(final_distance_max),' m' + write(IMAIN,*) ! add warning if estimate is poor ! (usually means receiver outside the mesh given by the user) if (final_distance_max > elemsize_max_glob) then - write(IMAIN,*) write(IMAIN,*) '************************************************************' write(IMAIN,*) '************************************************************' write(IMAIN,*) '***** WARNING: at least one receiver is poorly located *****' write(IMAIN,*) '************************************************************' write(IMAIN,*) '************************************************************' + write(IMAIN,*) endif ! write the locations of stations, so that we can load them and write them to SU headers later @@ -428,7 +436,6 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) ! elapsed time since beginning of mesh generation tCPU = wtime() - tstart - write(IMAIN,*) write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',tCPU write(IMAIN,*) write(IMAIN,*) 'End of receiver detection - done' diff --git a/src/specfem3D/read_stations.f90 b/src/specfem3D/read_stations.f90 index 2007c9d73..b12c9d578 100644 --- a/src/specfem3D/read_stations.f90 +++ b/src/specfem3D/read_stations.f90 @@ -38,12 +38,9 @@ subroutine read_stations(rec_filename) character(len=*),intent(in) :: rec_filename ! local parameters - integer :: i,irec,ier + integer :: irec,ier character(len=MAX_STRING_LEN) :: line - integer, allocatable, dimension(:) :: station_duplet - integer :: length_station_name,length_network_name - ! loop on all the stations to read the file if (myrank == 0) then @@ -67,20 +64,73 @@ subroutine read_stations(rec_filename) ! close receiver file close(IIN) - ! In case that the same station and network name appear twice (or more times) in the STATIONS - ! file, problems occur, as two (or more) seismograms are written (with mode - ! "append") to a file with same name. The philosophy here is to accept multiple - ! appearances and to just add a count to the station name in this case. - allocate(station_duplet(nrec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1970') - if (ier /= 0 ) call exit_MPI(myrank,'Error allocating station_duplet array') + ! find duplicate station names + call read_stations_find_duplets() + endif - station_duplet(:) = 0 - do irec = 1,nrec - ! checks if duplicate of another station read in so far + ! broadcast values to other slices + call bcast_all_ch_array(station_name,nrec,MAX_LENGTH_STATION_NAME) + call bcast_all_ch_array(network_name,nrec,MAX_LENGTH_NETWORK_NAME) + call bcast_all_dp(stlat,nrec) + call bcast_all_dp(stlon,nrec) + call bcast_all_dp(stele,nrec) + call bcast_all_dp(stbur,nrec) + + end subroutine read_stations + +! +!------------------------------------------------------------------------------------------- +! + + subroutine read_stations_find_duplets() + + use constants, only: MAX_LENGTH_STATION_NAME,MAX_LENGTH_NETWORK_NAME, & + DUPLETS_NREC_MINIMUM_FOR_HASH + + use specfem_par, only: myrank,nrec,station_name,network_name + + implicit none + + ! local parameters + integer :: i,irec,ier + + integer, allocatable, dimension(:) :: station_duplet + integer :: length_station_name,length_network_name + + ! hash table + integer, allocatable, dimension(:) :: hash_table + integer :: hash,hash_prob,hash_collisions + + ! In case that the same station and network name appear twice (or more times) in the STATIONS + ! file, problems occur, as two (or more) seismograms are written (with mode + ! "append") to a file with same name. The philosophy here is to accept multiple + ! appearances and to just add a count to the station name in this case. + allocate(station_duplet(nrec),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 1970') + if (ier /= 0 ) call exit_MPI(myrank,'Error allocating station_duplet array') + station_duplet(:) = 0 + + ! initializes the hash table + if (nrec >= DUPLETS_NREC_MINIMUM_FOR_HASH) then + ! makes hash table slightly larger than the actual number of stations to avoid many hash collisions + allocate(hash_table(5*nrec),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 1971') + if (ier /= 0 ) call exit_MPI(myrank,'Error allocating hash_table array') + hash_table(:) = 0 + hash_collisions = 0 + endif + + do irec = 1,nrec + ! checks if duplicate of another station read in so far + if (nrec < DUPLETS_NREC_MINIMUM_FOR_HASH) then + ! way 1: + ! loops over all previous stations and checks station/network names + ! this will get slow for very large STATIONS files (> 100,000 stations); scales approx. quadratic ~O(n**2) do i = 1,irec-1 if ((station_name(irec) == station_name(i)) .and. (network_name(irec) == network_name(i))) then + ! increases duplet count station_duplet(i) = station_duplet(i) + 1 + ! appends duplet number to station name if (len_trim(station_name(irec)) <= MAX_LENGTH_STATION_NAME-3) then write(station_name(irec),"(a,'_',i2.2)") trim(station_name(irec)),station_duplet(i)+1 else @@ -89,45 +139,153 @@ subroutine read_stations(rec_filename) endif enddo - ! checks name lengths - length_station_name = len_trim(station_name(irec)) - length_network_name = len_trim(network_name(irec)) + else + ! way 2: + ! gets a hash value and checks in hash table for duplicates; scales approx. linearly ~O(n) + hash = hashFunc(station_name(irec),network_name(irec)) + if (hash_table(hash) == 0) then + ! stores station index + hash_table(hash) = irec + else + ! found a duplicate hash + ! check if name matches + i = hash_table(hash) + if ((station_name(irec) == station_name(i)) .and. (network_name(irec) == network_name(i))) then + ! debug + !print *,'debug: Duplicate found in hash table:', & + ! irec,trim(station_name(irec)),trim(network_name(irec)), & + ! ' - hash number',hash, & + ! 'return index',i,trim(station_name(i)),trim(network_name(i)) + ! increases duplet count + station_duplet(i) = station_duplet(i) + 1 + ! appends duplet number to station name + if (len_trim(station_name(irec)) <= MAX_LENGTH_STATION_NAME-3) then + write(station_name(irec),"(a,'_',i2.2)") trim(station_name(irec)),station_duplet(i)+1 + else + call exit_MPI(myrank,'Please increase MAX_LENGTH_STATION_NAME by at least 3 to name station duplets') + endif + else + ! hash collision + hash_collisions = hash_collisions + 1 + ! debug + !print *,'debug: Collision found in hash table:', & + ! irec,trim(station_name(irec)),trim(network_name(irec)), & + ! ' - hash number',hash, & + ! 'return index',i,trim(station_name(i)),trim(network_name(i)) + ! put hash in next free slot (linear probing) + hash_prob = hash + do while (hash_table(hash_prob) /= 0) + ! increases hash index + hash_prob = mod(hash_prob + 1,size(hash_table)) + ! check if we reach again same hash, then table is full + if (hash_prob == hash) then + print *,'Error: Hash table is full, please consider a larger hash table!' + call exit_MPI(myrank,'Please increase hash table size for station duplets search') + endif + ! check entry + i = hash_table(hash_prob) + if (i == 0) then + ! stores station index in new free slot + hash_table(hash_prob) = irec + exit + else if (i == irec) then + ! station already stored, done adding + exit + else + ! check entry, maybe station moved hash index due to previous collisions + if ((station_name(irec) == station_name(i)) .and. (network_name(irec) == network_name(i))) then + ! debug + !print *,'debug: Duplicate found in hash table:', & + ! irec,trim(station_name(irec)),trim(network_name(irec)), & + ! ' - hash number by collision probing',hash_prob, & + ! 'return index',i,trim(station_name(i)),trim(network_name(i)) + ! increases duplet count + station_duplet(i) = station_duplet(i) + 1 + ! appends duplet number to station name + if (len_trim(station_name(irec)) <= MAX_LENGTH_STATION_NAME-3) then + write(station_name(irec),"(a,'_',i2.2)") trim(station_name(irec)),station_duplet(i)+1 + else + call exit_MPI(myrank,'Please increase MAX_LENGTH_STATION_NAME by at least 3 to name station duplets') + endif + ! done + exit + endif + endif + enddo + endif - ! check that length conforms to standard - if (length_station_name < 1 .or. length_station_name > MAX_LENGTH_STATION_NAME) then - print *, 'Error: invalid station name ',trim(station_name(irec)) - call exit_MPI(myrank,'wrong length of station name') - endif - if (length_network_name < 1 .or. length_network_name > MAX_LENGTH_NETWORK_NAME) then - print *, 'Error: invalid network name ',trim(network_name(irec)) - call exit_MPI(myrank,'wrong length of network name') endif + endif - enddo + ! checks name lengths + length_station_name = len_trim(station_name(irec)) + length_network_name = len_trim(network_name(irec)) - ! user output - if (maxval(station_duplet) > 0) then - print *,'Warning: found ',maxval(station_duplet),' station duplets (having same network & station names)' - print *,' station_duplet: ',station_duplet(:) - print *,' station_name : ',station_name(:) - print *,'Please check your STATIONS file entries to avoid confusions...' - print * + ! check that length conforms to standard + if (length_station_name < 1 .or. length_station_name > MAX_LENGTH_STATION_NAME) then + print *, 'Error: invalid station name ',trim(station_name(irec)) + call exit_MPI(myrank,'wrong length of station name') + endif + if (length_network_name < 1 .or. length_network_name > MAX_LENGTH_NETWORK_NAME) then + print *, 'Error: invalid network name ',trim(network_name(irec)) + call exit_MPI(myrank,'wrong length of network name') endif - ! free temporary array - deallocate(station_duplet) + enddo + ! user output + if (sum(station_duplet) > 0) then + print * + print *,'Warning: found ',sum(station_duplet),' station duplets (having same network & station names)' + do irec = 1,nrec + if (station_duplet(irec) > 0) then + print *,' station_name : ',station_name(irec),' network: ',network_name(irec) + print *,' irec : ',irec,' duplets: ',station_duplet(irec) + endif + enddo + print *,'Please check your STATIONS file entries to avoid confusions...' + print * endif - ! broadcast values to other slices - call bcast_all_ch_array(station_name,nrec,MAX_LENGTH_STATION_NAME) - call bcast_all_ch_array(network_name,nrec,MAX_LENGTH_NETWORK_NAME) - call bcast_all_dp(stlat,nrec) - call bcast_all_dp(stlon,nrec) - call bcast_all_dp(stele,nrec) - call bcast_all_dp(stbur,nrec) + ! debug: hash table info + !if (nrec >= DUPLETS_NREC_MINIMUM_FOR_HASH .and. hash_collisions > 0) & + ! print *,'debug: hash table collisions: ',hash_collisions + + ! free temporary array + deallocate(station_duplet) + if (nrec >= DUPLETS_NREC_MINIMUM_FOR_HASH) deallocate(hash_table) + + contains + + ! defines a simple hash function + integer function hashFunc(sta_name,net_name) + character(len=MAX_LENGTH_STATION_NAME), intent(in) :: sta_name + character(len=MAX_LENGTH_NETWORK_NAME), intent(in) :: net_name + ! local parameters + integer :: i, sum, prime + integer, parameter :: base = 31 ! seems to be a good number (arbitrary choice for a prime) + ! smaller prime numbers lead to more collisions + character(len=(MAX_LENGTH_STATION_NAME+MAX_LENGTH_NETWORK_NAME)) :: name + + ! full name, e.g., S00001DB + name = sta_name // net_name + + sum = 0 + prime = 1 + do i = 1, len_trim(name) + ! too simple, will get lots of hash collisions + ! sum = sum + ichar(name(i:i)) + + ! adding a bit more complexity, based on polynomial rolling + sum = mod(sum + ichar(name(i:i)) * prime, size(hash_table)) + prime = mod(prime * base, size(hash_table)) + enddo - end subroutine read_stations + hashFunc = mod(sum, size(hash_table)) + end function hashFunc + + + end subroutine read_stations_find_duplets ! !-------------------------------------------------------------------------------------------