updates writing seismograms

src/decompose_mesh/rules.mk

@@ -154,6 +154,7 @@ xdecompose_mesh_mpi_OBJECTS = \
 	$O/module_database.dec.o \
 	$O/module_mesh.dec.o \
 	$O/module_partition.dec.o \
+	$O/part_decompose_mesh.dec_module.o \
 	$O/program_decompose_mesh_mpi.mpidec.o \
 	$(EMPTY_MACRO)
 

src/specfem3D/compute_seismograms.f90

@@ -634,7 +634,7 @@ subroutine compute_seismograms_moment_adjoint()
                                       hprime_xx,hprime_yy,hprime_zz)
 
       ! source time function value
-      time_source_dble = dble(NSTEP-it)*DT - t0 - tshift_src(irec)
+      time_source_dble = dble(NSTEP-it) * DT - t0 - tshift_src(irec)
 
       ! determines source time function value
       stf = get_stf_viscoelastic(time_source_dble,irec,NSTEP-it+1)

src/specfem3D/setup_sources_receivers.f90

@@ -1265,13 +1265,27 @@ subroutine setup_receivers_precompute_intp()
   !         nadj_rec_local - determines the number of adjoint sources, i.e., number of station locations (STATIONS_ADJOINT), which
   !                          act as sources to drive the adjoint wavefield
 
+  ! check if we need to save seismos
+  if (SIMULATION_TYPE == 3 .and. (.not. SAVE_SEISMOGRAMS_IN_ADJOINT_RUN)) then
+    do_save_seismograms = .false.
+  else
+    do_save_seismograms = .true.
+  endif
+
+  ! sets local receivers to zero if no seismogram needs to be saved
+  if (.not. do_save_seismograms) nrec_local = 0
+
   ! user output
   if (myrank == 0) then
     write(IMAIN,*) 'seismograms:'
-    if (WRITE_SEISMOGRAMS_BY_MAIN) then
-      write(IMAIN,*) '  seismograms written by main process only'
+    if (do_save_seismograms) then
+      if (WRITE_SEISMOGRAMS_BY_MAIN) then
+        write(IMAIN,*) '  seismograms written by main process only'
+      else
+        write(IMAIN,*) '  seismograms written by all processes'
+      endif
     else
-      write(IMAIN,*) '  seismograms written by all processes'
+      write(IMAIN,*) '  seismograms will not be saved'
     endif
     write(IMAIN,*)
     write(IMAIN,*) '  Total number of simulation steps (NSTEP)                       = ',NSTEP

src/specfem3D/specfem3D_par.F90

@@ -193,6 +193,7 @@ module specfem_par
 
   integer :: nlength_seismogram
   integer :: seismo_offset,seismo_current
+  logical :: do_save_seismograms
 
   ! for ASDF/SAC headers time
   integer :: yr_PDE,jda_PDE,ho_PDE,mi_PDE

src/specfem3D/write_output_ASCII_or_binary.f90

@@ -37,14 +37,17 @@ subroutine write_output_ASCII_or_binary(one_seismogram, &
 
   use constants
 
-  use specfem_par, only: myrank,USE_BINARY_FOR_SEISMOGRAMS,SAVE_ALL_SEISMOS_IN_ONE_FILE,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
-    seismo_offset,seismo_current,NTSTEP_BETWEEN_OUTPUT_SAMPLE
+  use specfem_par, only: myrank, &
+    USE_BINARY_FOR_SEISMOGRAMS,SAVE_ALL_SEISMOS_IN_ONE_FILE, &
+    nlength_seismogram, &
+    NTSTEP_BETWEEN_OUTPUT_SAMPLE, &
+    seismo_offset,seismo_current
 
   implicit none
 
   integer,intent(in) :: NSTEP,it,SIMULATION_TYPE
 
-  real(kind=CUSTOM_REAL), dimension(NDIM,NTSTEP_BETWEEN_OUTPUT_SEISMOS/NTSTEP_BETWEEN_OUTPUT_SAMPLE),intent(in) :: one_seismogram
+  real(kind=CUSTOM_REAL), dimension(NDIM,nlength_seismogram),intent(in) :: one_seismogram
 
   double precision,intent(in) :: t0,DT
 
@@ -56,7 +59,6 @@ subroutine write_output_ASCII_or_binary(one_seismogram, &
   integer :: isample,ier,it_current
   double precision :: time_t_db
   real(kind=CUSTOM_REAL) :: value,time_t
-
   real, dimension(1:seismo_current) :: tr
 
   ! opens seismogram file
@@ -112,11 +114,9 @@ subroutine write_output_ASCII_or_binary(one_seismogram, &
     ! time
     if (SIMULATION_TYPE == 1) then
       ! forward simulation
-      ! distinguish between single and double precision for reals
-      time_t_db = dble( (it_current-1) * NTSTEP_BETWEEN_OUTPUT_SAMPLE) * DT - t0
+      time_t_db = dble((it_current-1) * NTSTEP_BETWEEN_OUTPUT_SAMPLE) * DT - t0
     else if (SIMULATION_TYPE == 3) then
       ! adjoint simulation: backward/reconstructed wavefields
-      ! distinguish between single and double precision for reals
       ! note: compare time_t with time used for source term
       time_t_db = dble(NSTEP - it_current * NTSTEP_BETWEEN_OUTPUT_SAMPLE) * DT - t0
     endif

src/specfem3D/write_output_ASDF.f90

@@ -156,8 +156,8 @@ end subroutine write_output_ASDF
 !! \param nrec_store The number of receivers on the local processor
   subroutine init_asdf_data(nrec_store)
 
-  use constants, only: NDIM
-  use specfem_par, only: myrank,NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS
+  use constants, only: NDIM,myrank
+  use specfem_par, only: NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS
 
   use asdf_data, only: asdf_container
 
@@ -221,7 +221,9 @@ subroutine store_asdf_data(seismogram_tmp, irec_local, irec, chn, iorientation)
 
   use constants, only: CUSTOM_REAL,NDIM,myrank
 
-  use specfem_par, only: NSTEP,NTSTEP_BETWEEN_OUTPUT_SAMPLE,nlength_seismogram, &
+  use specfem_par, only: &
+    nlength_seismogram, &
+    seismo_current, &
     stlat,stlon,stele,stbur,station_name,network_name
 
   use asdf_data, only: asdf_container
@@ -236,10 +238,6 @@ subroutine store_asdf_data(seismogram_tmp, irec_local, irec, chn, iorientation)
   ! local Variables
   integer :: length_station_name, length_network_name
   integer :: ier, i, index_increment
-  integer :: seismo_current_used
-
-  ! actual seismogram length
-  seismo_current_used = ceiling(real(NSTEP) / NTSTEP_BETWEEN_OUTPUT_SAMPLE)
 
   index_increment = iorientation
 
@@ -258,13 +256,13 @@ subroutine store_asdf_data(seismogram_tmp, irec_local, irec, chn, iorientation)
   asdf_container%receiver_el(irec_local) = stele(irec)
   asdf_container%receiver_dpt(irec_local) = stbur(irec)
 
-  allocate(asdf_container%records(i)%record(seismo_current_used),stat=ier)
+  allocate(asdf_container%records(i)%record(seismo_current),stat=ier)
   if (ier /= 0) call exit_MPI_without_rank('error allocating array 2017')
   if (ier /= 0) call exit_MPI (myrank, 'Allocating ASDF container failed.')
   asdf_container%records(i)%record(:) = 0.0
 
   ! seismogram as real data
-  asdf_container%records(i)%record(1:seismo_current_used) = real(seismogram_tmp(iorientation, 1:seismo_current_used),kind=4)
+  asdf_container%records(i)%record(1:seismo_current) = real(seismogram_tmp(iorientation, 1:seismo_current),kind=4)
 
   end subroutine store_asdf_data
 
@@ -316,7 +314,10 @@ subroutine write_asdf()
   use iso_c_binding, only: C_NULL_CHAR !,c_ptr
 !  use iso_Fortran_env
 
-  use specfem_par, only: seismo_offset,DT,NSTEP,NTSTEP_BETWEEN_OUTPUT_SAMPLE,OUTPUT_FILES,WRITE_SEISMOGRAMS_BY_MAIN
+  use specfem_par, only: DT,NSTEP, &
+    seismo_offset,seismo_current, &
+    NTSTEP_BETWEEN_OUTPUT_SAMPLE, &
+    OUTPUT_FILES,WRITE_SEISMOGRAMS_BY_MAIN
 
   implicit none
 
@@ -327,7 +328,6 @@ subroutine write_asdf()
   integer :: num_stations
   integer :: stationxml_length
   integer :: nsamples  ! constant, as in SPECFEM
-  integer :: seismo_current_used
   double precision :: sampling_rate
   double precision :: startTime
   integer(kind=8) :: start_time
@@ -393,9 +393,6 @@ subroutine write_asdf()
   call world_duplicate(comm)
   call world_size(mysize)
 
-  ! actual seismogram length
-  seismo_current_used = ceiling(real(NSTEP) / NTSTEP_BETWEEN_OUTPUT_SAMPLE)
-
   num_stations = asdf_container%nrec_local
   sampling_rate = 1.0/(DT*NTSTEP_BETWEEN_OUTPUT_SAMPLE)
 
@@ -620,7 +617,7 @@ subroutine write_asdf()
   deallocate(displs)
   deallocate(rcounts)
 
-  allocate(one_seismogram(NDIM,seismo_current_used),stat=ier)
+  allocate(one_seismogram(NDIM,seismo_current),stat=ier)
   if (ier /= 0) call exit_MPI_without_rank('error allocating array 2032')
   if (ier /= 0) call exit_MPI(myrank,'error allocating one_seismogram array')
   one_seismogram(:,:) = 0.0
@@ -795,7 +792,7 @@ subroutine write_asdf()
               if (len_trim(component_names_gather(i+((j-1)*NDIM),k)) == 0) cycle
 
               !  write(*,*) j, l, l+i, size(asdf_container%records)
-              one_seismogram(i,:) = asdf_container%records(l+i)%record(1:seismo_current_used)
+              one_seismogram(i,:) = asdf_container%records(l+i)%record(1:seismo_current)
             enddo
           endif
 
@@ -810,14 +807,14 @@ subroutine write_asdf()
               !       we therefore skip components with an empty ("") entry.
               if (len_trim(component_names_gather(i+((j-1)*NDIM),k)) == 0) cycle
 
-              one_seismogram(i,:) = asdf_container%records(l+i)%record(1:seismo_current_used)
+              one_seismogram(i,:) = asdf_container%records(l+i)%record(1:seismo_current)
             enddo
             ! send (real) data
-            call send_r(one_seismogram,NDIM*seismo_current_used,receiver,itag)
+            call send_r(one_seismogram,NDIM*seismo_current,receiver,itag)
 
           else if (myrank == 0) then
             ! receive (real) data
-            call recv_r(one_seismogram,NDIM*seismo_current_used,sender,itag)
+            call recv_r(one_seismogram,NDIM*seismo_current,sender,itag)
 
           endif
         endif
@@ -852,7 +849,7 @@ subroutine write_asdf()
 
             ! writes (float) data
             call ASDF_write_partial_waveform_f(data_ids(i), &
-                                               one_seismogram(i,1:seismo_current_used), 0, seismo_current_used, ier)
+                                               one_seismogram(i,1:seismo_current), 0, seismo_current, ier)
             if (ier /= 0) call exit_MPI(myrank,'Error ASDF write partial waveform failed')
 
             call ASDF_close_dataset_f(data_ids(i), ier)
@@ -938,11 +935,11 @@ subroutine write_asdf()
           ! see asdf implementation of ASDF_write_partial_waveform_f():
           ! https://github.com/SeismicData/asdf-library/blob/f28233894ef0a065eb725b9d22e55d0d02a7aac1/src/ASDF_write.c#L354
           if (myrank == current_proc) then
-            one_seismogram(i,:) = asdf_container%records(l+i)%record(1:seismo_current_used)
+            one_seismogram(i,:) = asdf_container%records(l+i)%record(1:seismo_current)
 
             ! writes (float) data
             call ASDF_write_partial_waveform_f(data_ids(i), &
-                                               one_seismogram(i,1:seismo_current_used), 0, seismo_current_used, ier)
+                                               one_seismogram(i,1:seismo_current), 0, seismo_current, ier)
             if (ier /= 0) call exit_MPI(myrank,'Error ASDF parallel write partial waveform failed')
           endif
 

src/specfem3D/write_seismograms.f90

@@ -40,14 +40,9 @@ subroutine write_seismograms()
   ! timing
   double precision, external :: wtime
   double precision :: write_time_begin,write_time
-  logical :: do_save_seismograms
 
-  ! check if we need to save seismos
-  if (SIMULATION_TYPE == 3 .and. (.not. SAVE_SEISMOGRAMS_IN_ADJOINT_RUN)) then
-    do_save_seismograms = .false.
-  else
-    do_save_seismograms = .true.
-  endif
+  ! checks if anything to do
+  if (.not. do_save_seismograms) return
 
   ! checks subsampling recurrence
   if (mod(it-1,NTSTEP_BETWEEN_OUTPUT_SAMPLE) == 0) then
@@ -75,7 +70,7 @@ subroutine write_seismograms()
     ! - ispec_selected_source   -> element containing source position, which becomes an adjoint "receiver"
 
     ! gets seismogram values
-    if (do_save_seismograms .and. nrec_local > 0) then
+    if (nrec_local > 0) then
       ! seismograms displ/veloc/accel/pressure
       if (GPU_MODE) then
         ! on GPU
@@ -98,9 +93,23 @@ subroutine write_seismograms()
   endif ! NTSTEP_BETWEEN_OUTPUT_SAMPLE
 
   ! write the current or final seismograms
+  !
+  ! for example:
+  !   NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10 , NTSTEP_BETWEEN_OUTPUT_SAMPLE = 5
+  !       -> 10 / 5 == 2 steps (nlength_seismogram)
+  !          we will store samples at it==1,it==6 and output at it==10
+  !   NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10 , NTSTEP_BETWEEN_OUTPUT_SAMPLE = 1
+  !       -> 10 / 1 == 10 steps
+  !          we will store samples at it==1,2,3,..,10 and output at it==10
+  !
+  !   that is for NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10 , NTSTEP_BETWEEN_OUTPUT_SAMPLE = 5,
+  !   the seismograms have samples for the wavefield at mod((it-1),NTSTEP_BETWEEN_OUTPUT_SAMPLE) == 0,
+  !   which is at it==1,it==6. for writing the seismograms however,
+  !   we would write out at mod(it,NTSTEP_BETWEEN_OUTPUT_SEISMOS) == 0, which is at it==10
+  !
   if (mod(it,NTSTEP_BETWEEN_OUTPUT_SEISMOS) == 0 .or. it == it_end) then
 
-    if (do_save_seismograms .and. .not. INVERSE_FWI_FULL_PROBLEM) then
+    if (.not. INVERSE_FWI_FULL_PROBLEM) then
 
       ! user output
       if (myrank == 0) then
@@ -157,7 +166,7 @@ subroutine write_seismograms()
         ! timing
         write_time = wtime() - write_time_begin
         ! output
-        write(IMAIN,*) 'Total number of time steps done: ', it-it_begin+1
+        write(IMAIN,*) 'Total number of time steps written: ', seismo_offset + seismo_current
         if (WRITE_SEISMOGRAMS_BY_MAIN) then
           write(IMAIN,*) 'Writing the seismograms by main proc alone took ',sngl(write_time),' seconds'
         else
@@ -167,7 +176,7 @@ subroutine write_seismograms()
         call flush_IMAIN()
       endif
 
-    endif ! do_save_seismograms
+    endif
 
     ! resets current seismogram position
     seismo_offset = seismo_offset + seismo_current
@@ -856,7 +865,7 @@ subroutine write_adj_seismograms_to_file(istore)
           it_tmp = seismo_offset + isample
 
           ! subtract half duration of the source to make sure travel time is correct
-          time_t = real(dble((it_tmp-1)*NTSTEP_BETWEEN_OUTPUT_SAMPLE) * DT - t0,kind=CUSTOM_REAL)
+          time_t = real(dble((it_tmp-1) * NTSTEP_BETWEEN_OUTPUT_SAMPLE) * DT - t0,kind=CUSTOM_REAL)
 
           ! distinguish between single and double precision for reals
           write(IOUT,*) time_t,all_seismograms(iorientation,isample,irec_local)