updates adjoint source time indexing

Author: danielpeter

src/specfem3D/compute_add_sources_acoustic.f90

@@ -211,6 +211,9 @@ subroutine compute_add_sources_acoustic(potential_dot_dot_acoustic)
 
       ! adds source term
       if (it < NSTEP) then
+        ! time step index for adjoint source (time-reversed)
+        it_index = NTSTEP_BETWEEN_READ_ADJSRC - mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC)
+
         ! receivers act as sources
         do irec_local = 1, nadj_rec_local
           irec = number_adjsources_global(irec_local)
@@ -225,9 +228,6 @@ subroutine compute_add_sources_acoustic(potential_dot_dot_acoustic)
 
                   hlagrange = hxir_adjstore(i,irec_local) * hetar_adjstore(j,irec_local) * hgammar_adjstore(k,irec_local)
 
-                  ! time step index
-                  it_index = NTSTEP_BETWEEN_READ_ADJSRC - mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC)
-
                   ! beware, for acoustic medium, a pressure source would be taking the negative
                   ! and divide by Kappa of the fluid;
                   ! this would have to be done when constructing the adjoint source.

src/specfem3D/compute_add_sources_poroelastic.f90

@@ -59,7 +59,7 @@ subroutine compute_add_sources_poroelastic()
   double precision,external :: get_stf_poroelastic
 
   logical :: ibool_read_adj_arrays
-  integer :: isource,iglob,i,j,k,ispec,it_sub_adj
+  integer :: isource,iglob,i,j,k,ispec,it_sub_adj,it_index
   integer :: irec_local,irec
   real(kind=CUSTOM_REAL) :: phil,tortl,rhol_s,rhol_f,rhol_bar
   real(kind=CUSTOM_REAL) :: fac_s,fac_w
@@ -247,6 +247,9 @@ subroutine compute_add_sources_poroelastic()
       endif ! if (ibool_read_adj_arrays)
 
       if (it < NSTEP) then
+        ! time step index for adjoint source (time-reversed)
+        it_index = NTSTEP_BETWEEN_READ_ADJSRC - mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC)
+
         ! receivers act as sources
         do irec_local = 1, nadj_rec_local
           irec = number_adjsources_global(irec_local)
@@ -271,12 +274,11 @@ subroutine compute_add_sources_poroelastic()
 
                   ! solid phase
                   accels_poroelastic(:,iglob) = accels_poroelastic(:,iglob) &
-                    + source_adjoint(:,irec_local,NTSTEP_BETWEEN_READ_ADJSRC - mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC)) * hlagrange
+                                              + source_adjoint(:,irec_local,it_index) * hlagrange
                   !
                   ! fluid phase
                   accelw_poroelastic(:,iglob) = accelw_poroelastic(:,iglob) &
-                    - rhol_f/rhol_bar &
-                    * source_adjoint(:,irec_local,NTSTEP_BETWEEN_READ_ADJSRC - mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC)) * hlagrange
+                                              - rhol_f/rhol_bar * source_adjoint(:,irec_local,it_index) * hlagrange
                 enddo
               enddo
             enddo
@@ -363,10 +365,10 @@ subroutine compute_add_sources_poroelastic()
 
                 ! solid phase
                 b_accels_poroelastic(:,iglob) = b_accels_poroelastic(:,iglob) &
-                             + fac_s * sourcearrays(isource,:,i,j,k)*stf_used
+                                              + fac_s * sourcearrays(isource,:,i,j,k)*stf_used
                 ! fluid phase
                 b_accelw_poroelastic(:,iglob) = b_accelw_poroelastic(:,iglob) &
-                             + fac_w * sourcearrays(isource,:,i,j,k)*stf_used
+                                              + fac_w * sourcearrays(isource,:,i,j,k)*stf_used
               enddo
             enddo
           enddo

src/specfem3D/compute_add_sources_viscoelastic.F90

@@ -70,7 +70,7 @@ subroutine compute_add_sources_viscoelastic(accel)
   double precision :: stf,time_source_dble,time_t
   double precision,external :: get_stf_viscoelastic
 
-  integer :: isource,iglob,i,j,k,ispec,it_sub_adj
+  integer :: isource,iglob,i,j,k,ispec,it_sub_adj,it_index
   integer :: irec_local,irec
 
   character(len=MAX_STRING_LEN) :: adj_source_file
@@ -230,6 +230,9 @@ subroutine compute_add_sources_viscoelastic(accel)
 
       ! adds source term
       if (it < NSTEP) then
+        ! time step index for adjoint source (time-reversed)
+        it_index = NTSTEP_BETWEEN_READ_ADJSRC - mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC)
+
         ! receivers act as sources
         do irec_local = 1, nadj_rec_local
           irec = number_adjsources_global(irec_local)
@@ -244,8 +247,7 @@ subroutine compute_add_sources_viscoelastic(accel)
 
                   hlagrange = hxir_adjstore(i,irec_local) * hetar_adjstore(j,irec_local) * hgammar_adjstore(k,irec_local)
 
-                  accel(:,iglob) = accel(:,iglob) &
-                    + source_adjoint(:,irec_local,NTSTEP_BETWEEN_READ_ADJSRC - mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC)) * hlagrange
+                  accel(:,iglob) = accel(:,iglob) + source_adjoint(:,irec_local,it_index) * hlagrange
                 enddo
               enddo
             enddo

src/specfem3D/compute_arrays_source.f90

@@ -213,7 +213,7 @@ subroutine compute_arrays_adjoint_source(adj_source_file,irec_local)
   character(len=*), intent(in) :: adj_source_file
 
   ! local
-  integer :: icomp, itime, ier, it_start, it_end, it_sub_adj
+  integer :: icomp, itime, ier, it_start, it_end, it_sub_adj, index_i
   real(kind=CUSTOM_REAL), dimension(NTSTEP_BETWEEN_READ_ADJSRC) :: adj_source_asdf
   real(kind=CUSTOM_REAL) :: val,junk
   ! note: should have same order as orientation in write_seismograms_to_file()
@@ -244,8 +244,12 @@ subroutine compute_arrays_adjoint_source(adj_source_file,irec_local)
 
       ! store the block we need
       do itime = it_start, it_end
+
+        ! index will run from 1 to NTSTEP_BETWEEN_READ_ADJSRC
+        index_i = itime - it_start + 1
+
         ! store adjoint trace
-        source_adjoint(icomp,irec_local,itime-it_start+1) = adj_source_asdf(itime-it_start+1)
+        source_adjoint(icomp,irec_local,index_i) = adj_source_asdf(index_i)
       enddo
     enddo
 
@@ -278,8 +282,11 @@ subroutine compute_arrays_adjoint_source(adj_source_file,irec_local)
           call exit_MPI(myrank,'file '//trim(filename)//' has wrong length, please check with your simulation duration (2)')
         endif
 
+        ! index will run from 1 to NTSTEP_BETWEEN_READ_ADJSRC
+        index_i = itime - it_start + 1
+
         ! store adjoint trace
-        source_adjoint(icomp,irec_local,itime-it_start+1) = val
+        source_adjoint(icomp,irec_local,index_i) = val
       enddo
 
       close(IIN)