Merge pull request #1673 from danielpeter/devel

updates reading/writing stations infos for large number of stations (>1000); updates and adds python scripts for setup and visualization

src/generate_databases/setup_mesh_adjacency.f90

@@ -220,8 +220,7 @@ subroutine setup_mesh_adjacency()
 
     ! user output
     if (myrank == 0) then
-      write(IMAIN,*) '     using kd-tree search radius = ',sngl(r_search)
-      write(IMAIN,*)
+      write(IMAIN,*) '     using kd-tree search radius             = ',sngl(r_search)
       call flush_IMAIN()
     endif
 
@@ -530,6 +529,17 @@ subroutine setup_mesh_adjacency()
       deallocate(kdtree_search_index)
     endif
 
+    ! user output progress
+    if (myrank == 0) then
+      if (mod(ispec_ref,max(NSPEC_AB/10,1)) == 0) then
+        tCPU = wtime() - time1
+        ! elapsed
+        write(IMAIN,*) "    ",int(ispec_ref/(max(NSPEC_AB/10,1)) * 10)," %", &
+                       " - elapsed time:",sngl(tCPU),"s"
+        ! flushes file buffer for main output file (IMAIN)
+        call flush_IMAIN()
+      endif
+    endif
   enddo ! ispec_ref
 
   ! frees temporary array
@@ -646,11 +656,12 @@ subroutine setup_mesh_adjacency()
   if (myrank == 0) then
     ! elapsed time since beginning of neighbor detection
     tCPU = wtime() - time1
+    write(IMAIN,*)
     write(IMAIN,*) '     maximum search elements                                      = ',num_elements_max
     write(IMAIN,*) '     maximum of actual search elements (after distance criterion) = ',num_elements_actual_max
     write(IMAIN,*)
     write(IMAIN,*) '     estimated maximum element size            = ',elemsize_max_glob
-    write(IMAIN,*) '     maximum distance between neighbor centers = ',sqrt(dist_squared_max)
+    write(IMAIN,*) '     maximum distance between neighbor centers = ',real(sqrt(dist_squared_max),kind=CUSTOM_REAL)
     if (use_kdtree_search) then
       if (sqrt(dist_squared_max) > r_search - 0.5*elemsize_max_glob) then
         write(IMAIN,*) '***'
@@ -665,7 +676,7 @@ subroutine setup_mesh_adjacency()
     endif
     write(IMAIN,*) '     total number of neighbors           = ',num_neighbors_all
     write(IMAIN,*)
-    write(IMAIN,*) '     Elapsed time for detection of neighbors in seconds = ',tCPU
+    write(IMAIN,*) '     Elapsed time for detection of neighbors in seconds = ',sngl(tCPU)
     write(IMAIN,*)
     call flush_IMAIN()
   endif

src/meshfem3D/check_mesh_quality.f90

@@ -348,7 +348,7 @@ subroutine check_mesh_quality(VP_MAX,NGLOB,NSPEC,x,y,z,ibool, &
     endif
 
     ! vtk file output
-    call write_VTK_data_elem_cr_meshfem(nspec,nglob,x,y,z,ibool,tmp1,filename)
+    call write_VTK_data_elem_cr_meshfem(nspec,nglob,NGLLX_M,x,y,z,ibool,tmp1,filename)
 
     deallocate(tmp1)
   endif

src/shared/write_VTK_data.f90

@@ -963,19 +963,19 @@ end subroutine write_VTK_data_ngnod_elem_i
 !
 !------------------------------------------------------------------------------------
 
-  subroutine write_VTK_data_elem_cr_meshfem(nspec,nglob,xstore_db,ystore_db,zstore_db,ibool, &
+  subroutine write_VTK_data_elem_cr_meshfem(nspec,nglob,NGLL,xstore_db,ystore_db,zstore_db,ibool, &
                                             elem_data,filename)
 
 ! special routine for meshfem3D with simpler mesh arrays
 
-  use constants, only: CUSTOM_REAL,MAX_STRING_LEN,NGNOD_EIGHT_CORNERS,IOUT_VTK
+  use constants, only: CUSTOM_REAL,MAX_STRING_LEN,IOUT_VTK
 
   implicit none
 
-  integer :: nspec,nglob
+  integer, intent(in) :: nspec,nglob,NGLL
 
   ! global coordinates
-  integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC) :: ibool
+  integer,dimension(NGLL,NGLL,NGLL,nspec),intent(in) :: ibool
   double precision, dimension(nglob) :: xstore_db,ystore_db,zstore_db
 
   ! element flag array
@@ -1004,8 +1004,10 @@ subroutine write_VTK_data_elem_cr_meshfem(nspec,nglob,xstore_db,ystore_db,zstore
   write(IOUT_VTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
   do ispec = 1,nspec
     write(IOUT_VTK,'(9i12)') 8, &
-                             ibool(1,ispec)-1,ibool(2,ispec)-1,ibool(4,ispec)-1,ibool(3,ispec)-1, &
-                             ibool(5,ispec)-1,ibool(6,ispec)-1,ibool(8,ispec)-1,ibool(7,ispec)-1
+                             ibool(1,1,1,ispec)-1,ibool(NGLL,1,1,ispec)-1, &
+                             ibool(NGLL,NGLL,1,ispec)-1,ibool(1,NGLL,1,ispec)-1, &
+                             ibool(1,1,NGLL,ispec)-1,ibool(NGLL,1,NGLL,ispec)-1, &
+                             ibool(NGLL,NGLL,NGLL,ispec)-1,ibool(1,NGLL,NGLL,ispec)-1
   enddo
   write(IOUT_VTK,*) ''
 

src/specfem3D/locate_receivers.f90

@@ -32,7 +32,7 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source)
 
   use constants
 
-  use specfem_par, only: USE_SOURCES_RECEIVERS_Z,SUPPRESS_UTM_PROJECTION,INVERSE_FWI_FULL_PROBLEM,SU_FORMAT, &
+  use specfem_par, only: USE_SOURCES_RECEIVERS_Z,SUPPRESS_UTM_PROJECTION,INVERSE_FWI_FULL_PROBLEM, &
                          ibool,myrank,NSPEC_AB,NGLOB_AB, &
                          xstore,ystore,zstore, &
                          nrec,islice_selected_rec,ispec_selected_rec, &
@@ -113,6 +113,7 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source)
     if (USE_SOURCES_RECEIVERS_Z) then
       write(IMAIN,*) 'using sources/receivers Z:'
       write(IMAIN,*) '  (depth) becomes directly (z) coordinate'
+      write(IMAIN,*)
     endif
     call flush_IMAIN()
 
@@ -132,16 +133,18 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source)
   call read_stations(rec_filename)
 
   ! checks if station locations already available
-  if (SU_FORMAT .and. (.not. INVERSE_FWI_FULL_PROBLEM) ) then
-    call read_stations_SU_from_previous_run(is_done_stations)
+  if (.not. INVERSE_FWI_FULL_PROBLEM) then
+    ! reads stations_info.bin if available in OUTPUT_FILES/ folder
+    call read_stations_from_previous_run(is_done_stations)
+
     ! check if done
     if (is_done_stations) then
       ! user output
       if (myrank == 0) then
         ! elapsed time since beginning of mesh generation
         tCPU = wtime() - tstart
         write(IMAIN,*)
-        write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',tCPU
+        write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',sngl(tCPU)
         write(IMAIN,*)
         write(IMAIN,*) 'End of receiver detection - done'
         write(IMAIN,*)
@@ -317,7 +320,7 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source)
       endif
 
       ! limits user output if too many receivers
-      if (nrec < 1000 .and. (.not. SU_FORMAT )) then
+      if (nrec < 1000) then
 
         ! receiver info
         write(IMAIN,*)
@@ -432,11 +435,11 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source)
     close(IOUT_SU)
 
     ! stores station infos for later runs
-    if (SU_FORMAT) call write_stations_SU_for_next_run(x_found,y_found,z_found)
+    call write_stations_for_next_run(x_found,y_found,z_found)
 
     ! elapsed time since beginning of mesh generation
     tCPU = wtime() - tstart
-    write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',tCPU
+    write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',sngl(tCPU)
     write(IMAIN,*)
     write(IMAIN,*) 'End of receiver detection - done'
     write(IMAIN,*)

src/specfem3D/locate_source.F90

@@ -703,7 +703,7 @@ subroutine locate_source()
   if (myrank == 0) then
     tCPU = wtime() - tstart
     write(IMAIN,*)
-    write(IMAIN,*) 'Elapsed time for detection of sources in seconds = ',tCPU
+    write(IMAIN,*) 'Elapsed time for detection of sources in seconds = ',sngl(tCPU)
     write(IMAIN,*)
     write(IMAIN,*) 'End of source detection - done'
     write(IMAIN,*)