Skip to content

Commit

Permalink
updates adjoint source reading (ASDF & SU format)
Browse files Browse the repository at this point in the history
  • Loading branch information
@danielpeter
danielpeter committed Nov 13, 2024
1 parent a0ae26d commit b4d8af3
Showing 1 changed file with 133 additions and 62 deletions.
195 changes: 133 additions & 62 deletions src/specfem3D/compute_arrays_source.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -214,9 +214,8 @@ subroutine compute_arrays_adjoint_source(adj_source_file,irec_local)

! local
integer :: icomp, itime, ier, it_start, it_end, it_sub_adj
real(kind=CUSTOM_REAL), dimension(NDIM,NTSTEP_BETWEEN_READ_ADJSRC) :: adj_src
real(kind=CUSTOM_REAL), dimension(NSTEP) :: adj_source_asdf
double precision :: junk
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()
character(len=3),dimension(NDIM) :: comp
character(len=MAX_STRING_LEN) :: filename
Expand All @@ -231,31 +230,30 @@ subroutine compute_arrays_adjoint_source(adj_source_file,irec_local)
it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1

adj_src(:,:) = 0._CUSTOM_REAL
itime = 0

if (READ_ADJSRC_ASDF) then
! ASDF format
do icomp = 1,NDIM ! 3 components
! format: "net_sta_comp"
filename = trim(adj_source_file) // '_' // comp(icomp)

! would skip read and set source artificially to zero if out of bounds,
! see comments above
if (it_start == 0 .and. itime == 0) then
adj_src(icomp,1) = 0._CUSTOM_REAL
cycle
endif

! reads full trace (NSTEP)
call read_adjoint_sources_ASDF(filename, adj_source_asdf, it_start, it_end)

! stores source array
adj_src(icomp,:) = adj_source_asdf(:)
! debug - check whether we read the correct block
!if (icomp == 1) print *, junk, adj_source_asdf(itime-it_start+1,icomp)

! store the block we need
do itime = it_start, it_end
! store adjoint trace
source_adjoint(icomp,irec_local,itime-it_start+1) = adj_source_asdf(itime-it_start+1)
enddo
enddo

else
! ASCII format
! loops over components
do icomp = 1, NDIM
! format: "SEM/net.sta.comp.adj"
filename = OUTPUT_FILES(1:len_trim(OUTPUT_FILES))//'/../SEM/'//trim(adj_source_file)//'.'//comp(icomp)//'.adj'

open(unit=IIN,file=trim(filename),status='old',action='read',iostat = ier)
Expand All @@ -268,19 +266,20 @@ subroutine compute_arrays_adjoint_source(adj_source_file,irec_local)
!! skip unused blocks
do itime = 1, it_start-1
read(IIN,*,iostat=ier) junk, junk
if (ier /= 0) &
call exit_MPI(myrank, &
'file '//trim(filename)//' has wrong length, please check with your simulation duration (1111)')
if (ier /= 0) then
call exit_MPI(myrank,'file '//trim(filename)//' has wrong length, please check with your simulation duration (1)')
endif
enddo

!! read the block we need
do itime = it_start, it_end
read(IIN,*,iostat=ier) junk, source_adjoint(icomp,irec_local,itime-it_start+1)
!!! used to check whether we read the correct block
! if (icomp==1) print *, junk, adj_src(itime-it_start+1,icomp)
if (ier /= 0) &
call exit_MPI(myrank, &
'file '//trim(filename)//' has wrong length, please check with your simulation duration (2222)')
read(IIN,*,iostat=ier) junk, val
if (ier /= 0) then
call exit_MPI(myrank,'file '//trim(filename)//' has wrong length, please check with your simulation duration (2)')
endif

! store adjoint trace
source_adjoint(icomp,irec_local,itime-it_start+1) = val
enddo

close(IIN)
Expand All @@ -294,72 +293,144 @@ end subroutine compute_arrays_adjoint_source
!-------------------------------------------------------------------------------------------------
!

subroutine compute_arrays_adjoint_source_SU()
subroutine compute_arrays_adjoint_source_SU(idomain)

use specfem_par, only: myrank,source_adjoint,it,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC,nrec_local
use shared_parameters, only: ACOUSTIC_SIMULATION,ELASTIC_SIMULATION
use constants

implicit none
integer, intent(in) :: idomain
! local parameters
real(kind=CUSTOM_REAL), dimension(NTSTEP_BETWEEN_READ_ADJSRC) :: adj_temp
integer :: ier, irec_local, it_start, it_sub_adj
logical :: found_adjoint_files

! note: should have same order as orientation in write_seismograms_to_file()
character(len=MAX_STRING_LEN) :: procname, filename
character(len=MAX_STRING_LEN) :: procname, filename_p, filename_x, filename_y, filename_z

! check if anything to read for this slice
if (nrec_local < 1) return

! range of the block we need to read
it_sub_adj = ceiling( dble(it)/dble(NTSTEP_BETWEEN_READ_ADJSRC) )
it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC

write(procname,"(i4)") myrank

if (ACOUSTIC_SIMULATION) then
! check if we have adjoint traces
found_adjoint_files = .false.

select case(idomain)
case (IDOMAIN_ACOUSTIC)
! get acoustic adjoint traces
! SU adjoint file name
filename_p = trim(OUTPUT_FILES)//'../SEM/'//trim(adjustl(procname))//'_p_SU.adj'

! check if file exists
inquire(file=trim(filename_p),exist=found_adjoint_files)

! read file
if (found_adjoint_files) then
! user output
if (myrank == 0) then
write(IMAIN,*) 'reading acoustic adjoint traces:'
write(IMAIN,*) ' ',trim(filename_p)
write(IMAIN,*) ' using SU_FORMAT'
write(IMAIN,*)
write(IMAIN,*) ' start index = ',it_start
write(IMAIN,*) ' trace length = ',NTSTEP_BETWEEN_READ_ADJSRC
write(IMAIN,*)
call flush_IMAIN()
endif

filename = trim(OUTPUT_FILES)//'../SEM/'//trim(adjustl(procname))//'_p_SU.adj'
open(unit=IIN_SU1,file=filename,status='old',access='stream',iostat = ier)
if (ier /= 0) call exit_MPI(myrank,'file '//trim(filename)//' does not exist')
do irec_local = 1,nrec_local
read(IIN_SU1,pos=4*((irec_local-1)*(60+NSTEP) + 60 + it_start)+1 ) adj_temp
source_adjoint(1,irec_local,:) = adj_temp(:)
source_adjoint(2,irec_local,:) = 0.0 !TRIVIAL
source_adjoint(3,irec_local,:) = 0.0 !TRIVIAL
enddo
close(IIN_SU1)
open(unit=IIN_SU1,file=trim(filename_p),status='old',access='stream',iostat = ier)
if (ier /= 0) call exit_MPI(myrank,'file '//trim(filename_p)//' does not exist')
do irec_local = 1,nrec_local
read(IIN_SU1,pos=4*((irec_local-1)*(60+NSTEP) + 60 + it_start)+1 ) adj_temp
source_adjoint(1,irec_local,:) = adj_temp(:)
source_adjoint(2,irec_local,:) = 0.0_CUSTOM_REAL !TRIVIAL
source_adjoint(3,irec_local,:) = 0.0_CUSTOM_REAL !TRIVIAL
enddo
close(IIN_SU1)
endif

case (IDOMAIN_ELASTIC)
! get elastic adjoint traces
! SU adjoint file names
! x-direction traces
filename_x = trim(OUTPUT_FILES)//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj'
! y-direction traces
filename_y = trim(OUTPUT_FILES)//'../SEM/'//trim(adjustl(procname))//'_dy_SU.adj'
! z-direction traces
filename_z = trim(OUTPUT_FILES)//'../SEM/'//trim(adjustl(procname))//'_dz_SU.adj'

! check if x-file exists
inquire(file=trim(filename_x),exist=found_adjoint_files)
if (found_adjoint_files) then
! check if y-file exists
inquire(file=trim(filename_y),exist=found_adjoint_files)
if (found_adjoint_files) then
! check if z-file exists
inquire(file=trim(filename_z),exist=found_adjoint_files)
endif
endif

if (found_adjoint_files) then
! user output
if (myrank == 0) then
write(IMAIN,*) 'reading elastic adjoint traces:'
write(IMAIN,*) ' ',trim(filename_x)
write(IMAIN,*) ' ',trim(filename_y)
write(IMAIN,*) ' ',trim(filename_z)
write(IMAIN,*) ' using SU_FORMAT'
write(IMAIN,*)
write(IMAIN,*) ' start index = ',it_start
write(IMAIN,*) ' trace length = ',NTSTEP_BETWEEN_READ_ADJSRC
write(IMAIN,*)
call flush_IMAIN()
endif

else if (ELASTIC_SIMULATION) then
! opens files
open(unit=IIN_SU1,file=trim(filename_x),status='old',access='stream',iostat = ier)
if (ier /= 0) call exit_MPI(myrank,'file '//trim(filename_x)//' does not exist')

filename = trim(OUTPUT_FILES)//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj'
open(unit=IIN_SU1,file=filename,status='old',access='stream',iostat = ier)
if (ier /= 0) call exit_MPI(myrank,'file '//trim(filename)//' does not exist')
open(unit=IIN_SU2,file=trim(filename_y),status='old',access='stream',iostat = ier)
if (ier /= 0) call exit_MPI(myrank,'file '//trim(filename_y)//' does not exist')

filename = trim(OUTPUT_FILES)//'../SEM/'//trim(adjustl(procname))//'_dy_SU.adj'
open(unit=IIN_SU2,file=filename,status='old',access='stream',iostat = ier)
if (ier /= 0) call exit_MPI(myrank,'file '//trim(filename)//' does not exist')
open(unit=IIN_SU3,file=trim(filename_z),status='old',access='stream',iostat = ier)
if (ier /= 0) call exit_MPI(myrank,'file '//trim(filename_z)//' does not exist')

filename = trim(OUTPUT_FILES)//'../SEM/'//trim(adjustl(procname))//'_dz_SU.adj'
open(unit=IIN_SU3,file=filename,status='old',access='stream',iostat = ier)
if (ier /= 0) call exit_MPI(myrank,'file '//trim(filename)//' does not exist')
! reads traces
do irec_local = 1,nrec_local
read(IIN_SU1,pos=4*((irec_local-1)*(60+NSTEP) + 60 + it_start)+1 ) adj_temp
source_adjoint(1,irec_local,:) = adj_temp(:)

do irec_local = 1,nrec_local
read(IIN_SU1,pos=4*((irec_local-1)*(60+NSTEP) + 60 + it_start)+1 ) adj_temp
source_adjoint(1,irec_local,:) = adj_temp(:)
read(IIN_SU2,pos=4*((irec_local-1)*(60+NSTEP) + 60 + it_start)+1 ) adj_temp
source_adjoint(2,irec_local,:) = adj_temp(:)

read(IIN_SU2,pos=4*((irec_local-1)*(60+NSTEP) + 60 + it_start)+1 ) adj_temp
source_adjoint(2,irec_local,:) = adj_temp(:)
read(IIN_SU3,pos=4*((irec_local-1)*(60+NSTEP) + 60 + it_start)+1 ) adj_temp
source_adjoint(3,irec_local,:) = adj_temp(:)
enddo

read(IIN_SU3,pos=4*((irec_local-1)*(60+NSTEP) + 60 + it_start)+1 ) adj_temp
source_adjoint(3,irec_local,:) = adj_temp(:)
enddo
close(IIN_SU1)
close(IIN_SU2)
close(IIN_SU3)
endif

close(IIN_SU1)
close(IIN_SU2)
close(IIN_SU3)
case (IDOMAIN_POROELASTIC)
! not implemented yet
call exit_MPI(myrank,'SU_FORMAT not implemented for adjoint sources in poroelastic domains yet')

else
case default
! domain not recognized
call exit_MPI(myrank,'Invalid domain for SU_FORMAT adjoint sources')

call exit_MPI(myrank,'SU_FORMAT not implemented for adjoint poroelastic simulations yet')
end select

endif
! debug - check if file found
!if (.not. found_adjoint_files) then
! call exit_MPI(myrank,'Found no adjoint traces in SU_FORMAT')
!endif

end subroutine compute_arrays_adjoint_source_SU

Expand Down

0 comments on commit b4d8af3

Please sign in to comment.