Skip to content

Commit

Permalink
updates adding source contributions
Browse files Browse the repository at this point in the history
  • Loading branch information
@danielpeter
danielpeter committed Nov 13, 2024
1 parent b4d8af3 commit 16db5a1
Show file tree
Hide file tree
Showing 6 changed files with 342 additions and 288 deletions.
117 changes: 69 additions & 48 deletions src/specfem3D/compute_add_sources_acoustic.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,9 @@ subroutine compute_add_sources_acoustic(potential_dot_dot_acoustic)

use constants
use specfem_par, only: station_name,network_name, &
nsources_local,tshift_src,DT,t0,SU_FORMAT,USE_LDDRK,istage, &
nsources_local,tshift_src,DT,t0, &
SU_FORMAT,READ_ADJSRC_ASDF, &
USE_LDDRK,istage, &
hxir_adjstore,hetar_adjstore,hgammar_adjstore,source_adjoint,number_adjsources_global,nadj_rec_local, &
USE_BINARY_FOR_SEISMOGRAMS, &
ibool,NSOURCES,myrank,it,ispec_selected_source,islice_selected_source, &
Expand Down Expand Up @@ -62,27 +64,28 @@ subroutine compute_add_sources_acoustic(potential_dot_dot_acoustic)

character(len=MAX_STRING_LEN) :: adj_source_file

! sets current initial time
if (USE_LDDRK) then
! LDDRK
! note: the LDDRK scheme updates displacement after the stiffness computations and
! 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(it-1-1)*DT + dble(C_LDDRK(istage))*DT - t0
else
time_t = dble(it-1)*DT - t0
endif

! forward simulations
! forward simulations
if (SIMULATION_TYPE == 1 .and. nsources_local > 0) then

! ignore pressure sources for fault rupture simulations
if (FAULT_SIMULATION) return

! no source inside the mesh if we are coupling with DSM
! because the source is precisely the wavefield coming from the DSM traction file
if (COUPLE_WITH_INJECTION_TECHNIQUE) return

! sets current initial time
if (USE_LDDRK) then
! LDDRK
! note: the LDDRK scheme updates displacement after the stiffness computations and
! 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(it-1-1)*DT + dble(C_LDDRK(istage))*DT - t0
else
time_t = dble(it-1)*DT - t0
endif

! openmp solver
!$OMP PARALLEL if (NSOURCES > 100) &
!$OMP DEFAULT(SHARED) &
Expand Down Expand Up @@ -179,25 +182,34 @@ subroutine compute_add_sources_acoustic(potential_dot_dot_acoustic)
! with other partitions while we calculate for the inner part
! this must be done carefully, otherwise the adjoint sources may be added twice
if (ibool_read_adj_arrays .and. .not. INVERSE_FWI_FULL_PROBLEM) then

if (.not. SU_FORMAT) then
! ASCII format
! reads adjoint source files
if (SU_FORMAT) then
! SU format
call compute_arrays_adjoint_source_SU(IDOMAIN_ACOUSTIC)
else if (READ_ADJSRC_ASDF) then
! ASDF format
do irec_local = 1, nadj_rec_local
! reads in **net**.**sta**.**BH**.adj files
irec = number_adjsources_global(irec_local)
adj_source_file = trim(network_name(irec))//'_'//trim(station_name(irec)) ! format: "net_sta"
! compute source arrays
call compute_arrays_adjoint_source(adj_source_file,irec_local)
enddo
else
! default ASCII format
if (USE_BINARY_FOR_SEISMOGRAMS) stop 'Adjoint simulations not supported with .bin format, please use SU format instead'

!!! read ascii adjoint sources
do irec_local = 1, nadj_rec_local
irec = number_adjsources_global(irec_local)
! reads in **net**.**sta**.**BH**.adj files
adj_source_file = trim(network_name(irec))//'.'//trim(station_name(irec))
call compute_arrays_adjoint_source(adj_source_file,irec)
irec = number_adjsources_global(irec_local)
adj_source_file = trim(network_name(irec))//'.'//trim(station_name(irec)) ! format: "net.sta"
! compute source arrays
call compute_arrays_adjoint_source(adj_source_file,irec_local)
enddo
else
! SU format
call compute_arrays_adjoint_source_SU()
endif !if (.not. SU_FORMAT)

endif
endif ! if (ibool_read_adj_arrays)

! adds source term
if (it < NSTEP) then
! receivers act as sources
do irec_local = 1, nadj_rec_local
Expand Down Expand Up @@ -251,7 +263,8 @@ end subroutine compute_add_sources_acoustic
subroutine compute_add_sources_acoustic_backward(b_potential_dot_dot_acoustic)

use constants
use specfem_par, only: nsources_local,tshift_src,DT,t0,USE_LDDRK,istage, &
use specfem_par, only: nsources_local,tshift_src,DT,t0, &
USE_LDDRK,istage, &
ibool,NSOURCES,myrank,it,islice_selected_source,ispec_selected_source, &
sourcearrays,kappastore,SIMULATION_TYPE,NSTEP,NGLOB_AB

Expand Down Expand Up @@ -408,13 +421,14 @@ subroutine compute_add_sources_acoustic_GPU()

use constants
use specfem_par, only: station_name,network_name, &
nsources_local,tshift_src,DT,t0,SU_FORMAT,USE_LDDRK,istage, &
nsources_local,tshift_src,DT,t0, &
SU_FORMAT,READ_ADJSRC_ASDF, &
USE_LDDRK,istage, &
source_adjoint,nadj_rec_local,number_adjsources_global, &
USE_BINARY_FOR_SEISMOGRAMS, &
NSOURCES,it,SIMULATION_TYPE,NSTEP,nrec, &
NTSTEP_BETWEEN_READ_ADJSRC,Mesh_pointer, &
INVERSE_FWI_FULL_PROBLEM,run_number_of_the_source, &
GPU_MODE
INVERSE_FWI_FULL_PROBLEM,run_number_of_the_source

! coupling
use shared_parameters, only: COUPLE_WITH_INJECTION_TECHNIQUE
Expand All @@ -438,11 +452,9 @@ subroutine compute_add_sources_acoustic_GPU()

character(len=MAX_STRING_LEN) :: adj_source_file

! checks if anything to do
if (.not. GPU_MODE) return

! forward simulations
if (SIMULATION_TYPE == 1 .and. nsources_local > 0) then

! ignore pressure sources for fault rupture simulations
if (FAULT_SIMULATION) return

Expand Down Expand Up @@ -525,24 +537,34 @@ subroutine compute_add_sources_acoustic_GPU()
! with other partitions while we calculate for the inner part
! this must be done carefully, otherwise the adjoint sources may be added twice
if (ibool_read_adj_arrays .and. .not. INVERSE_FWI_FULL_PROBLEM) then

if (.not. SU_FORMAT) then
! ASCII format
! reads adjoint source files
if (SU_FORMAT) then
! SU format
call compute_arrays_adjoint_source_SU(IDOMAIN_ACOUSTIC)
else if (READ_ADJSRC_ASDF) then
! ASDF format
do irec_local = 1, nadj_rec_local
! reads in **net**.**sta**.**BH**.adj files
irec = number_adjsources_global(irec_local)
adj_source_file = trim(network_name(irec))//'_'//trim(station_name(irec)) ! format: "net_sta"
! compute source arrays
call compute_arrays_adjoint_source(adj_source_file,irec_local)
enddo
else
! default ASCII format
if (USE_BINARY_FOR_SEISMOGRAMS) stop 'Adjoint simulations not supported with .bin format, please use SU format instead'
!!! read ascii adjoint sources
do irec_local = 1, nadj_rec_local
irec = number_adjsources_global(irec_local)
! reads in **net**.**sta**.**BH**.adj files
adj_source_file = trim(network_name(irec))//'.'//trim(station_name(irec))
call compute_arrays_adjoint_source(adj_source_file,irec)
irec = number_adjsources_global(irec_local)
adj_source_file = trim(network_name(irec))//'.'//trim(station_name(irec)) ! format: "net.sta"
! compute source arrays
call compute_arrays_adjoint_source(adj_source_file,irec_local)
enddo
else
! SU format
call compute_arrays_adjoint_source_SU()
endif !if (.not. SU_FORMAT)

endif
endif ! if (ibool_read_adj_arrays)

! adds source term
if (it < NSTEP) then
! receivers act as sources
! on GPU
Expand All @@ -560,10 +582,11 @@ end subroutine compute_add_sources_acoustic_GPU
subroutine compute_add_sources_acoustic_backward_GPU()

use constants
use specfem_par, only: nsources_local,tshift_src,DT,t0,USE_LDDRK,istage, &
use specfem_par, only: nsources_local,tshift_src,DT,t0, &
USE_LDDRK,istage, &
NSOURCES,myrank,it, &
SIMULATION_TYPE,NSTEP, &
GPU_MODE,Mesh_pointer,run_number_of_the_source
Mesh_pointer,run_number_of_the_source
! undo_att
use specfem_par, only: UNDO_ATTENUATION_AND_OR_PML,NSUBSET_ITERATIONS,NT_DUMP_ATTENUATION, &
iteration_on_subset,it_of_this_subset
Expand All @@ -587,8 +610,6 @@ subroutine compute_add_sources_acoustic_backward_GPU()
! checks if anything to do
if (SIMULATION_TYPE /= 3) return

if (.not. GPU_MODE) return

! checks if this slice has sources to add
if (nsources_local == 0) return

Expand Down
38 changes: 28 additions & 10 deletions src/specfem3D/compute_add_sources_poroelastic.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ subroutine compute_add_sources_poroelastic()
UNDO_ATTENUATION_AND_OR_PML, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source, &
sourcearrays,SIMULATION_TYPE,NSTEP, &
SU_FORMAT,READ_ADJSRC_ASDF, &
ispec_selected_rec, &
nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC, &
hxir_adjstore,hetar_adjstore,hgammar_adjstore,source_adjoint,number_adjsources_global,nadj_rec_local
Expand Down Expand Up @@ -67,6 +68,7 @@ subroutine compute_add_sources_poroelastic()

! forward simulations
if (SIMULATION_TYPE == 1 .and. nsources_local > 0) then

! ignore CMT sources for fault rupture simulations
if (FAULT_SIMULATION) return

Expand Down Expand Up @@ -217,16 +219,31 @@ subroutine compute_add_sources_poroelastic()
! this must be done carefully, otherwise the adjoint sources may be added
! twice
if (ibool_read_adj_arrays) then
if (USE_BINARY_FOR_SEISMOGRAMS) stop 'Adjoint simulations not supported with .bin format, please use ASCII instead'
! ASCII format
!!! read ascii adjoint sources
do irec_local = 1, nadj_rec_local
irec = number_adjsources_global(irec_local)
! compute source arrays
! reads in **net**.**sta**.**BH**.adj files
adj_source_file = trim(network_name(irec))//'.'//trim(station_name(irec))
call compute_arrays_adjoint_source(adj_source_file,irec_local)
enddo
! reads adjoint source files
if (SU_FORMAT) then
! SU format
call compute_arrays_adjoint_source_SU(IDOMAIN_ELASTIC)
else if (READ_ADJSRC_ASDF) then
! ASDF format
do irec_local = 1, nadj_rec_local
! reads in **net**.**sta**.**BH**.adj files
irec = number_adjsources_global(irec_local)
adj_source_file = trim(network_name(irec))//'_'//trim(station_name(irec)) ! format: "net_sta"
! compute source arrays
call compute_arrays_adjoint_source(adj_source_file,irec_local)
enddo
else
! default ASCII format
if (USE_BINARY_FOR_SEISMOGRAMS) stop 'Adjoint simulations not supported with .bin format, please use SU format instead'
!!! read ascii adjoint sources
do irec_local = 1, nadj_rec_local
! reads in **net**.**sta**.**BH**.adj files
irec = number_adjsources_global(irec_local)
adj_source_file = trim(network_name(irec))//'.'//trim(station_name(irec)) ! format: "net.sta"
! compute source arrays
call compute_arrays_adjoint_source(adj_source_file,irec_local)
enddo
endif
endif ! if (ibool_read_adj_arrays)

if (it < NSTEP) then
Expand Down Expand Up @@ -277,6 +294,7 @@ subroutine compute_add_sources_poroelastic()

! adjoint/backward simulations
if (SIMULATION_TYPE == 3 .and. nsources_local > 0) then

! ignore CMT sources for fault rupture simulations
if (FAULT_SIMULATION) return

Expand Down
Loading

0 comments on commit 16db5a1

Please sign in to comment.