From b4d8af3b6021f00d2991fe177f41323136addc27 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Wed, 13 Nov 2024 11:58:44 +0100 Subject: [PATCH] updates adjoint source reading (ASDF & SU format) --- src/specfem3D/compute_arrays_source.f90 | 195 ++++++++++++++++-------- 1 file changed, 133 insertions(+), 62 deletions(-) diff --git a/src/specfem3D/compute_arrays_source.f90 b/src/specfem3D/compute_arrays_source.f90 index 3d835744a..e8bc4d3a3 100644 --- a/src/specfem3D/compute_arrays_source.f90 +++ b/src/specfem3D/compute_arrays_source.f90 @@ -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 @@ -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) @@ -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) @@ -294,18 +293,23 @@ 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) ) @@ -313,53 +317,120 @@ subroutine compute_arrays_adjoint_source_SU() 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