Merge pull request #1756 from danielpeter/devel

updates point search to prefer slices that contain point locations inside elements

src/generate_databases/model_external_values.F90

@@ -94,7 +94,7 @@ subroutine model_external_broadcast()
   if (myrank == 0) call read_external_model()
 
   ! broadcast the information read on the main to the nodes
-  call bcast_all_dp(MEXT_V%dvs, size(MEXT_V%dvs))
+  call bcast_all_dp(MEXT_V%dvs, size(MEXT_V%dvs,kind=4))
 
   end subroutine model_external_broadcast
 

src/generate_databases/model_salton_trough.f90

@@ -71,7 +71,7 @@ subroutine model_salton_trough_broadcast()
   if (myrank == 0) call read_salton_sea_model()
 
   ! broadcast the information read on the main to the nodes
-  call bcast_all_r(vp_array, size(vp_array))
+  call bcast_all_r(vp_array, size(vp_array,kind=4))
 
   end subroutine model_salton_trough_broadcast
 

src/generate_databases/model_tomography.f90

@@ -121,7 +121,7 @@ subroutine model_tomography_broadcast(myrank)
   ! allocate( vp_tomography(1:nrecord) ,stat=ier)
   ! if (ier /= 0) stop 'error allocating array vp_tomography'
   !endif
-  !call bcast_all_cr(vp_tomography,size(vp_tomography))
+  !call bcast_all_cr(vp_tomography,size(vp_tomography,kind=4))
 
   ! synchronizes processes
   call synchronize_all()

src/generate_databases/setup_mesh_adjacency.f90

@@ -31,7 +31,7 @@ subroutine setup_mesh_adjacency()
 ! setups mesh adjacency array to search element neighbors for point searches
 
   use constants, only: myrank, &
-    NDIM,NGLLX,NGLLY,NGLLZ,MIDX,MIDY,MIDZ,IMAIN,CUSTOM_REAL,MAX_STRING_LEN
+    NGLLX,NGLLY,NGLLZ,MIDX,MIDY,MIDZ,IMAIN,MAX_STRING_LEN
 
   use generate_databases_par, only: NSPEC_AB,NGLOB_AB,ibool,NPROC,prname
 
@@ -255,7 +255,8 @@ subroutine setup_mesh_adjacency()
     else
       ! no neighbors
       ! warning
-      print *,'*** Warning: found mesh element with no neighbors : slice ',myrank,' - element ',ispec_ref,' ***'
+      print *,'*** Warning: found mesh element with no neighbors : slice ',myrank, &
+              ' - element ',ispec_ref,'out of',NSPEC_AB,' ***'
     endif
 
     ! again loop to get neighbors of neighbors
@@ -377,7 +378,7 @@ subroutine setup_mesh_adjacency()
   endif
 
   ! check if element has neighbors
-  ! note: in case of a fault in this slice (splitting nodes) and/or scotch paritioning
+  ! note: in case of a fault in this slice (splitting nodes) and/or scotch partitioning
   !       it can happen that an element has no neighbors
   if (NPROC == 1 .and. (.not. ANY_FAULT_IN_THIS_PROC)) then
     ! checks if neighbors were found

src/meshfem3D/determine_cavity.f90

@@ -366,13 +366,13 @@ subroutine cmm_determine_cavity(nglob)
         allocate(tmp_all(1,1),stat=ier)
         if (ier /= 0) call exit_MPI_without_rank('error allocating array 1330')
       endif
-      call sum_all_1Darray_dp(cavity_boundary,tmp_all,size(cavity_boundary))
+      call sum_all_1Darray_dp(cavity_boundary,tmp_all,size(cavity_boundary,kind=4))
       if (myrank == 0) then
         cavity_boundary(:,:) = tmp_all(:,:)
       endif
       deallocate(tmp_all)
       ! broadcasts to all others
-      call bcast_all_dp(cavity_boundary,size(cavity_boundary))
+      call bcast_all_dp(cavity_boundary,size(cavity_boundary,kind=4))
 
       !print *,'cavity boundary after:',myrank,'array:',cavity_boundary(:,:)
 

src/shared/recompute_jacobian.f90

@@ -30,7 +30,7 @@
   subroutine recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
                                 xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,NGNOD)
 
-  use constants, only: NDIM,ZERO
+  use constants, only: NDIM,ZERO,myrank
 
   implicit none
 
@@ -56,7 +56,7 @@ subroutine recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
 
 ! recompute jacobian for any (xi,eta,gamma) point, not necessarily a GLL point
 
-! check that the parameter file is correct
+  ! check that the parameter file is correct
   if (NGNOD /= 8 .and. NGNOD /= 27) stop 'elements should have 8 or 27 control nodes'
 
   if (NGNOD == 8) then
@@ -71,51 +71,66 @@ subroutine recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
   x = ZERO
   y = ZERO
   z = ZERO
+
   xxi = ZERO
   xeta = ZERO
   xgamma = ZERO
+
   yxi = ZERO
   yeta = ZERO
   ygamma = ZERO
+
   zxi = ZERO
   zeta = ZERO
   zgamma = ZERO
 
   do ia = 1,NGNOD
-    x = x+shape3D(ia)*xelm(ia)
-    y = y+shape3D(ia)*yelm(ia)
-    z = z+shape3D(ia)*zelm(ia)
+    x = x + shape3D(ia) * xelm(ia)
+    y = y + shape3D(ia) * yelm(ia)
+    z = z + shape3D(ia) * zelm(ia)
 
     ! Jacobian matrix elements
     ! dx/dxi, dx/deta, etc.
-    xxi = xxi+dershape3D(1,ia)*xelm(ia)
-    xeta = xeta+dershape3D(2,ia)*xelm(ia)
-    xgamma = xgamma+dershape3D(3,ia)*xelm(ia)
-    yxi = yxi+dershape3D(1,ia)*yelm(ia)
-    yeta = yeta+dershape3D(2,ia)*yelm(ia)
-    ygamma = ygamma+dershape3D(3,ia)*yelm(ia)
-    zxi = zxi+dershape3D(1,ia)*zelm(ia)
-    zeta = zeta+dershape3D(2,ia)*zelm(ia)
-    zgamma = zgamma+dershape3D(3,ia)*zelm(ia)
+    xxi = xxi + dershape3D(1,ia)*xelm(ia)
+    xeta = xeta + dershape3D(2,ia)*xelm(ia)
+    xgamma = xgamma + dershape3D(3,ia)*xelm(ia)
+
+    yxi = yxi + dershape3D(1,ia)*yelm(ia)
+    yeta = yeta + dershape3D(2,ia)*yelm(ia)
+    ygamma = ygamma + dershape3D(3,ia)*yelm(ia)
+
+    zxi = zxi + dershape3D(1,ia)*zelm(ia)
+    zeta = zeta + dershape3D(2,ia)*zelm(ia)
+    zgamma = zgamma + dershape3D(3,ia)*zelm(ia)
   enddo
 
   ! Jacobian determinant
   jacobian = xxi*(yeta*zgamma-ygamma*zeta) - xeta*(yxi*zgamma-ygamma*zxi) + xgamma*(yxi*zeta-yeta*zxi)
 
-  if (jacobian <= ZERO) stop '3D Jacobian undefined'
+  if (jacobian <= ZERO) then
+    ! info
+    print *,'Error Jacobian rank:',myrank,'has invalid Jacobian ',jacobian
+    print *,'  input xi/eta/gamma: ',xi,eta,gamma
+    print *,'  Jacobian: ',jacobian,'xxi',xxi,xeta,xgamma,'yxi',yxi,yeta,ygamma,'zxi',zxi,zeta,zgamma
+    print *,'  location x/y/z: ',x,y,z
+
+    stop '3D Jacobian undefined'
+  endif
 
   ! inverse Jacobian matrix
   !
   ! invert the relation (Fletcher p. 50 vol. 2)
-  xix = (yeta*zgamma-ygamma*zeta)/jacobian
-  xiy = (xgamma*zeta-xeta*zgamma)/jacobian
-  xiz = (xeta*ygamma-xgamma*yeta)/jacobian
-  etax = (ygamma*zxi-yxi*zgamma)/jacobian
-  etay = (xxi*zgamma-xgamma*zxi)/jacobian
-  etaz = (xgamma*yxi-xxi*ygamma)/jacobian
-  gammax = (yxi*zeta-yeta*zxi)/jacobian
-  gammay = (xeta*zxi-xxi*zeta)/jacobian
-  gammaz = (xxi*yeta-xeta*yxi)/jacobian
+  xix = (yeta*zgamma - ygamma*zeta) / jacobian
+  xiy = (xgamma*zeta - xeta*zgamma) / jacobian
+  xiz = (xeta*ygamma - xgamma*yeta) / jacobian
+
+  etax = (ygamma*zxi - yxi*zgamma) / jacobian
+  etay = (xxi*zgamma - xgamma*zxi) / jacobian
+  etaz = (xgamma*yxi - xxi*ygamma) / jacobian
+
+  gammax = (yxi*zeta - yeta*zxi) / jacobian
+  gammay = (xeta*zxi - xxi*zeta) / jacobian
+  gammaz = (xxi*yeta - xeta*yxi) / jacobian
 
   end subroutine recompute_jacobian
 

src/specfem3D/iterate_time.F90

@@ -418,9 +418,9 @@ subroutine it_transfer_from_GPU()
 
       if (ATTENUATION) then
         call transfer_rmemory_from_device(Mesh_pointer,R_xx,R_yy,R_xy,R_xz,R_yz, &
-                                          R_trace,size(R_xx))
+                                          R_trace,size(R_xx,kind=4))
         call transfer_strain_from_device(Mesh_pointer,epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
-                                         epsilondev_trace,size(epsilondev_xx))
+                                         epsilondev_trace,size(epsilondev_xx,kind=4))
 
       endif
     endif

src/specfem3D/locate_MPI_slice.f90

@@ -69,6 +69,17 @@ subroutine locate_MPI_slice(npoints_subset,ipoin_already_done, &
   double precision, dimension(npoints_subset,0:NPROC-1) :: final_distance_all
   double precision, dimension(NDIM,NDIM,npoints_subset,0:NPROC-1) :: nu_all
 
+  double precision :: xi,eta,gamma
+
+  ! point inside element
+  logical :: found_point_inside
+
+  ! prefer slice with point location inside element rather than just selecting slice with closest location
+  logical, parameter :: DO_PREFER_INSIDE_ELEMENT = .true.
+
+  !debug
+  !integer,dimension(1) :: iproc_min
+
   ! initializes with dummy values
   ispec_selected_all(:,:) = -1
   idomain_all(:,:) = -1000
@@ -98,33 +109,117 @@ subroutine locate_MPI_slice(npoints_subset,ipoin_already_done, &
   ! find the slice and element to put the source
   if (myrank == 0) then
 
-    ! loops over subset
-    do ipoin_in_this_subset = 1,npoints_subset
+    ! we prefer slices with point locations inside an element over slices with a closer point distance
+    ! but a point location outside the (best) element
+    if (DO_PREFER_INSIDE_ELEMENT) then
+
+      ! loops over subset
+      do ipoin_in_this_subset = 1,npoints_subset
+        ! mapping from station/source number in current subset to real station/source number in all the subsets
+        ipoin = ipoin_in_this_subset + ipoin_already_done
+
+        ! checks first if we have a close point inside an element
+        found_point_inside = .false.
+        distmin = HUGEVAL
+        do iproc = 0,NPROC-1
+          ! point position
+          xi = xi_all(ipoin_in_this_subset,iproc)
+          eta = eta_all(ipoin_in_this_subset,iproc)
+          gamma = gamma_all(ipoin_in_this_subset,iproc)
+
+          ! points inside an element have xi/eta/gamma <= 1.0
+          if (abs(xi) <= 1.d0 .and. abs(eta) <= 1.d0 .and. abs(gamma) <= 1.d0) then
+            ! takes point if closer
+            if (final_distance_all(ipoin_in_this_subset,iproc) < distmin) then
+              found_point_inside = .true.
+
+              xi_point(ipoin)    = xi
+              eta_point(ipoin)   = eta
+              gamma_point(ipoin) = gamma
+
+              distmin =  final_distance_all(ipoin_in_this_subset,iproc)
+
+              islice_selected(ipoin) = iproc
+              ispec_selected(ipoin) = ispec_selected_all(ipoin_in_this_subset,iproc)
+              idomain(ipoin) = idomain_all(ipoin_in_this_subset,iproc)
+
+              nu_point(:,:,ipoin) = nu_all(:,:,ipoin_in_this_subset,iproc)
+
+              x_found(ipoin) = x_found_all(ipoin_in_this_subset,iproc)
+              y_found(ipoin) = y_found_all(ipoin_in_this_subset,iproc)
+              z_found(ipoin) = z_found_all(ipoin_in_this_subset,iproc)
+            endif
+          endif
+        enddo
+
+        ! if we haven't found a close point inside an element, then look for just the closest possible
+        if (.not. found_point_inside) then
+          do iproc = 0,NPROC-1
+            if (final_distance_all(ipoin_in_this_subset,iproc) < distmin) then
+              distmin =  final_distance_all(ipoin_in_this_subset,iproc)
+
+              islice_selected(ipoin) = iproc
+              ispec_selected(ipoin) = ispec_selected_all(ipoin_in_this_subset,iproc)
+              idomain(ipoin) = idomain_all(ipoin_in_this_subset,iproc)
+
+              xi_point(ipoin)    = xi_all(ipoin_in_this_subset,iproc)
+              eta_point(ipoin)   = eta_all(ipoin_in_this_subset,iproc)
+              gamma_point(ipoin) = gamma_all(ipoin_in_this_subset,iproc)
+
+              nu_point(:,:,ipoin) = nu_all(:,:,ipoin_in_this_subset,iproc)
+
+              x_found(ipoin) = x_found_all(ipoin_in_this_subset,iproc)
+              y_found(ipoin) = y_found_all(ipoin_in_this_subset,iproc)
+              z_found(ipoin) = z_found_all(ipoin_in_this_subset,iproc)
+            endif
+          enddo
+        endif
 
-      ! mapping from station/source number in current subset to real station/source number in all the subsets
-      ipoin = ipoin_in_this_subset + ipoin_already_done
+        final_distance(ipoin) = distmin
+
+        !debug
+        !iproc_min = minloc(final_distance_all(ipoin_in_this_subset,:)) - 1 ! -1 to start procs at 0
+        !print *,'debug: locate MPI point',ipoin,ipoin_in_this_subset,'dist',final_distance_all(ipoin_in_this_subset,:), &
+        !  'iproc min',iproc_min(1), &
+        !  'xi min',xi_all(ipoin_in_this_subset,iproc_min(1)), &
+        !           eta_all(ipoin_in_this_subset,iproc_min(1)), &
+        !           gamma_all(ipoin_in_this_subset,iproc_min(1)), &
+        !  'iproc found',islice_selected(ipoin),'dist found',final_distance(ipoin), &
+        !  'xi found',xi_point(ipoin),eta_point(ipoin),gamma_point(ipoin)
+      enddo
 
-      distmin = HUGEVAL
-      do iproc = 0,NPROC-1
-        if (final_distance_all(ipoin_in_this_subset,iproc) < distmin) then
-          distmin =  final_distance_all(ipoin_in_this_subset,iproc)
+    else
+      ! old version takes closest point
 
-          islice_selected(ipoin) = iproc
-          ispec_selected(ipoin) = ispec_selected_all(ipoin_in_this_subset,iproc)
-          idomain(ipoin) = idomain_all(ipoin_in_this_subset,iproc)
+      ! loops over subset
+      do ipoin_in_this_subset = 1,npoints_subset
 
-          xi_point(ipoin)    = xi_all(ipoin_in_this_subset,iproc)
-          eta_point(ipoin)   = eta_all(ipoin_in_this_subset,iproc)
-          gamma_point(ipoin) = gamma_all(ipoin_in_this_subset,iproc)
-          nu_point(:,:,ipoin) = nu_all(:,:,ipoin_in_this_subset,iproc)
+        ! mapping from station/source number in current subset to real station/source number in all the subsets
+        ipoin = ipoin_in_this_subset + ipoin_already_done
 
-          x_found(ipoin) = x_found_all(ipoin_in_this_subset,iproc)
-          y_found(ipoin) = y_found_all(ipoin_in_this_subset,iproc)
-          z_found(ipoin) = z_found_all(ipoin_in_this_subset,iproc)
-        endif
+        distmin = HUGEVAL
+        do iproc = 0,NPROC-1
+          if (final_distance_all(ipoin_in_this_subset,iproc) < distmin) then
+            distmin =  final_distance_all(ipoin_in_this_subset,iproc)
+
+            islice_selected(ipoin) = iproc
+            ispec_selected(ipoin) = ispec_selected_all(ipoin_in_this_subset,iproc)
+            idomain(ipoin) = idomain_all(ipoin_in_this_subset,iproc)
+
+            xi_point(ipoin)    = xi_all(ipoin_in_this_subset,iproc)
+            eta_point(ipoin)   = eta_all(ipoin_in_this_subset,iproc)
+            gamma_point(ipoin) = gamma_all(ipoin_in_this_subset,iproc)
+            nu_point(:,:,ipoin) = nu_all(:,:,ipoin_in_this_subset,iproc)
+
+            x_found(ipoin) = x_found_all(ipoin_in_this_subset,iproc)
+            y_found(ipoin) = y_found_all(ipoin_in_this_subset,iproc)
+            z_found(ipoin) = z_found_all(ipoin_in_this_subset,iproc)
+          endif
+        enddo
+        final_distance(ipoin) = distmin
       enddo
-      final_distance(ipoin) = distmin
-    enddo
+
+    endif
 
   endif ! end of section executed by main process only