updates undo_attenuation timing for LDDRK

src/shared/get_attenuation_model.f90

@@ -34,10 +34,10 @@ module attenuation_model
     sequence
     double precision, dimension(:,:), pointer :: tau_eps_storage
     double precision, dimension(:), pointer :: Q_storage
-    integer Q_resolution
-    integer Q_max
+    integer :: Q_resolution
+    integer :: Q_max
   end type model_attenuation_storage_var
-  type (model_attenuation_storage_var) AM_S
+  type (model_attenuation_storage_var) :: AM_S
 
   ! attenuation_simplex_variables
   type attenuation_simplex_variables
@@ -756,7 +756,7 @@ subroutine get_attenuation_source_freq(f_c_source,min_period,max_period)
   f2 = 1.0d0 / min_period
 
   ! use the logarithmic central frequency
-  f_c_source = 10.0d0**(0.5 * (log10(f1) + log10(f2)))
+  f_c_source = 10.0d0**(0.5d0 * (log10(f1) + log10(f2)))
 
   end subroutine get_attenuation_source_freq
 

src/specfem3D/assemble_MPI_vector.f90

@@ -213,16 +213,14 @@ subroutine synchronize_MPI_vector_blocking_ord(NPROC,NGLOB_AB,array_val, &
     enddo
 
     ! set the value to that of the highest-rank processor
-  do iinterface = 1, num_interfaces_ext_mesh
-    if (myrank < my_neighbors_ext_mesh(iinterface)) then
-      do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
-        iglob = ibool_interfaces_ext_mesh(ipoin,iinterface)
-        array_val(:,iglob) = buffer_recv_vector(:,ipoin,iinterface)
-      enddo
-    endif
-  enddo
-
-
+    do iinterface = 1, num_interfaces_ext_mesh
+      if (myrank < my_neighbors_ext_mesh(iinterface)) then
+        do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+          iglob = ibool_interfaces_ext_mesh(ipoin,iinterface)
+          array_val(:,iglob) = buffer_recv_vector(:,ipoin,iinterface)
+        enddo
+      endif
+    enddo
 
   ! wait for communications completion (send)
     do iinterface = 1, num_interfaces_ext_mesh

src/specfem3D/compute_add_sources_acoustic.f90

@@ -317,7 +317,12 @@ subroutine compute_add_sources_acoustic_backward()
     !       after adding boundary/coupling/source terms.
     !       thus, at each time loop step it, displ(:) is still at (n) and not (n+1) like for the Newmark scheme
     !       when entering this routine. we therefore at an additional -DT to have the corresponding timing for the source.
-    time_t = dble(NSTEP-it_tmp-1)*DT - dble(C_LDDRK(istage))*DT - t0
+    if (UNDO_ATTENUATION_AND_OR_PML) then
+      ! stepping moves forward from snapshot position
+      time_t = dble(NSTEP-it_tmp-1)*DT + dble(C_LDDRK(istage))*DT - t0
+    else
+      time_t = dble(NSTEP-it_tmp-1)*DT - dble(C_LDDRK(istage))*DT - t0
+    endif
   else
     time_t = dble(NSTEP-it_tmp)*DT - t0
   endif
@@ -606,7 +611,13 @@ subroutine compute_add_sources_acoustic_backward_GPU()
     !       after adding boundary/coupling/source terms.
     !       thus, at each time loop step it, displ(:) is still at (n) and not (n+1) like for the Newmark scheme
     !       when entering this routine. we therefore at an additional -DT to have the corresponding timing for the source.
-    time_t = dble(NSTEP-it_tmp-1)*DT - dble(C_LDDRK(istage))*DT - t0
+    if (UNDO_ATTENUATION_AND_OR_PML) then
+      ! steps forward in time
+      time_t = dble(NSTEP-it_tmp-1)*DT + dble(C_LDDRK(istage))*DT - t0
+    else
+      ! steps backward
+      time_t = dble(NSTEP-it_tmp-1)*DT - dble(C_LDDRK(istage))*DT - t0
+    endif
   else
     time_t = dble(NSTEP-it_tmp)*DT - t0
   endif

src/specfem3D/compute_add_sources_poroelastic.f90

@@ -33,6 +33,7 @@ subroutine compute_add_sources_poroelastic()
   use specfem_par, only: station_name,network_name,USE_FORCE_POINT_SOURCE, &
                          tshift_src,dt,t0,USE_LDDRK,istage,USE_EXTERNAL_SOURCE_FILE,user_source_time_function, &
                          USE_BINARY_FOR_SEISMOGRAMS,ibool, &
+                         UNDO_ATTENUATION_AND_OR_PML, &
                          NSOURCES,myrank,it,islice_selected_source,ispec_selected_source, &
                          sourcearrays,SIMULATION_TYPE,NSTEP, &
                          ispec_selected_rec, &
@@ -287,7 +288,12 @@ subroutine compute_add_sources_poroelastic()
             !       after adding boundary/coupling/source terms.
             !       thus, at each time loop step it, displ(:) is still at (n) and not (n+1) like for the Newmark scheme
             !       when entering this routine. we therefore at an additional -DT to have the corresponding timing for the source.
-            time_source_dble = dble(NSTEP-it-1)*DT - dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+            if (UNDO_ATTENUATION_AND_OR_PML) then
+              ! stepping moves forward from snapshot position
+              time_source_dble = dble(NSTEP-it-1)*DT + dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+            else
+              time_source_dble = dble(NSTEP-it-1)*DT - dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+            endif
           else
             time_source_dble = dble(NSTEP-it)*DT - t0 - tshift_src(isource)
           endif

src/specfem3D/compute_add_sources_viscoelastic.F90

@@ -361,7 +361,12 @@ subroutine compute_add_sources_viscoelastic_backward()
     !       after adding boundary/coupling/source terms.
     !       thus, at each time loop step it, displ(:) is still at (n) and not (n+1) like for the Newmark scheme
     !       when entering this routine. we therefore at an additional -DT to have the corresponding timing for the source.
-    time_t = dble(NSTEP-it_tmp-1)*DT - dble(C_LDDRK(istage))*DT - t0
+    if (UNDO_ATTENUATION_AND_OR_PML) then
+      ! stepping moves forward from snapshot position
+      time_t = dble(NSTEP-it_tmp-1)*DT + dble(C_LDDRK(istage))*DT - t0
+    else
+      time_t = dble(NSTEP-it_tmp-1)*DT - dble(C_LDDRK(istage))*DT - t0
+    endif
   else
     time_t = dble(NSTEP-it_tmp)*DT - t0
   endif
@@ -442,6 +447,7 @@ subroutine compute_add_sources_viscoelastic_GPU()
                         num_free_surface_faces, &
                         nsources_local,tshift_src,dt,t0,SU_FORMAT, &
                         USE_LDDRK,istage,USE_EXTERNAL_SOURCE_FILE,user_source_time_function,USE_BINARY_FOR_SEISMOGRAMS, &
+                        UNDO_ATTENUATION_AND_OR_PML, &
                         NSOURCES,it,SIMULATION_TYPE,NSTEP, &
                         nrec,islice_selected_rec, &
                         nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
@@ -603,7 +609,12 @@ subroutine compute_add_sources_viscoelastic_GPU()
           !       after adding boundary/coupling/source terms.
           !       thus, at each time loop step it, displ(:) is still at (n) and not (n+1) like for the Newmark scheme
           !       when entering this routine. we therefore at an additional -DT to have the corresponding timing for the source.
-          time_source_dble = dble(NSTEP-it-1)*DT - dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+          if (UNDO_ATTENUATION_AND_OR_PML) then
+            ! stepping moves forward from snapshot position
+            time_source_dble = dble(NSTEP-it-1)*DT + dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+          else
+            time_source_dble = dble(NSTEP-it-1)*DT - dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+          endif
         else
           time_source_dble = dble(NSTEP-it)*DT - t0 - tshift_src(isource)
         endif
@@ -724,7 +735,12 @@ subroutine compute_add_sources_viscoelastic_backward_GPU()
       !       after adding boundary/coupling/source terms.
       !       thus, at each time loop step it, displ(:) is still at (n) and not (n+1) like for the Newmark scheme
       !       when entering this routine. we therefore at an additional -DT to have the corresponding timing for the source.
-      time_t = dble(NSTEP-it_tmp-1)*DT - dble(C_LDDRK(istage))*DT - t0
+      if (UNDO_ATTENUATION_AND_OR_PML) then
+        ! stepping moves forward from snapshot position
+        time_t = dble(NSTEP-it_tmp-1)*DT + dble(C_LDDRK(istage))*DT - t0
+      else
+        time_t = dble(NSTEP-it_tmp-1)*DT - dble(C_LDDRK(istage))*DT - t0
+      endif
     else
       time_t = dble(NSTEP-it_tmp)*DT - t0
     endif

src/specfem3D/prepare_attenuation.f90

@@ -98,8 +98,10 @@ subroutine prepare_attenuation()
 
   ! broadcasts
   call bcast_all_i_for_database(ispec, 1)
-  if (size(factor_common) > 0) call bcast_all_cr_for_database(factor_common(1,1,1,1,1), size(factor_common))
-  if (size(scale_factor) > 0) call bcast_all_cr_for_database(scale_factor(1,1,1,1), size(scale_factor))
+  if (size(factor_common) > 0) &
+    call bcast_all_cr_for_database(factor_common(1,1,1,1,1), size(factor_common))
+  if (size(scale_factor) > 0) &
+    call bcast_all_cr_for_database(scale_factor(1,1,1,1), size(scale_factor))
   call bcast_all_cr_for_database(factor_common_kappa(1,1,1,1,1), size(factor_common_kappa))
   call bcast_all_cr_for_database(scale_factor_kappa(1,1,1,1), size(scale_factor_kappa))
 

src/specfem3D/print_stf_file.f90

@@ -33,7 +33,7 @@ subroutine print_stf_file()
 
   use specfem_par, only: NSTEP,SIMULATION_TYPE,OUTPUT_FILES, &
     NSOURCES,islice_selected_source,ispec_selected_source,tshift_src,t0, &
-    DT,USE_LDDRK, &
+    DT,USE_LDDRK,UNDO_ATTENUATION_AND_OR_PML, &
     USE_EXTERNAL_SOURCE_FILE,user_source_time_function
 
   use specfem_par_acoustic, only: ispec_is_acoustic
@@ -98,7 +98,12 @@ subroutine print_stf_file()
             !       thus, at each time loop step it, displ(:) is still at (n) and not (n+1) like for the Newmark scheme
             !       when entering this routine. we therefore at an additional -DT to have the corresponding timing for the source.
             istage = 1  ! only 1. stage output (corresponds to u(n))
-            time_source_dble = dble(NSTEP-it-1)*DT - dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+            if (UNDO_ATTENUATION_AND_OR_PML) then
+              ! stepping moves forward from snapshot position
+              time_source_dble = dble(NSTEP-it-1)*DT + dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+            else
+              time_source_dble = dble(NSTEP-it-1)*DT - dble(C_LDDRK(istage))*DT - t0 - tshift_src(isource)
+            endif
           else
             time_source_dble = dble(NSTEP-it)*DT - t0 - tshift_src(isource)
           endif