[cig-commits] r22194 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Sun Jun 9 08:22:30 PDT 2013
Author: xie.zhinan
Date: 2013-06-09 08:22:30 -0700 (Sun, 09 Jun 2013)
New Revision: 22194
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_viscoelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_viscoelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90
Log:
finish splitting the part of code for simulation of backward wave propagation in stage of compute_forces_acoustic and compute_forces_viscoelastic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2013-06-08 12:12:16 UTC (rev 22193)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2013-06-09 15:22:30 UTC (rev 22194)
@@ -32,10 +32,9 @@
xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
- SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ SIMULATION_TYPE,NSTEP, &
nrec,islice_selected_rec,ispec_selected_rec, &
- nadj_rec_local,adj_sourcearrays,b_potential_dot_dot_acoustic, &
- NTSTEP_BETWEEN_READ_ADJSRC)
+ nadj_rec_local,adj_sourcearrays,NTSTEP_BETWEEN_READ_ADJSRC)
use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
@@ -74,11 +73,10 @@
logical, dimension(NSPEC_AB) :: ispec_is_acoustic
!adjoint simulations
- integer:: SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT
+ integer:: SIMULATION_TYPE,NSTEP
integer:: nrec
integer,dimension(nrec) :: islice_selected_rec,ispec_selected_rec
integer:: nadj_rec_local
- real(kind=CUSTOM_REAL),dimension(NGLOB_ADJOINT):: b_potential_dot_dot_acoustic
logical :: ibool_read_adj_arrays
integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: &
@@ -359,6 +357,104 @@
endif ! nadj_rec_local > 0
endif
+ ! master prints out source time function to file
+ if(PRINT_SOURCE_TIME_FUNCTION .and. phase_is_inner) then
+ time_source = (it-1)*DT - t0
+ call sum_all_cr(stf_used_total,stf_used_total_all)
+ if( myrank == 0 ) write(IOSTF,*) time_source,stf_used_total_all
+ endif
+
+ end subroutine compute_add_sources_acoustic
+
+!
+!=====================================================================
+! for acoustic solver for back propagation wave field
+
+ subroutine compute_add_sources_acoustic_bpwf(NSPEC_AB, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
+ xi_source,eta_source,gamma_source, &
+ hdur,hdur_gaussian,tshift_src,dt,t0, &
+ sourcearrays,kappastore,ispec_is_acoustic,&
+ SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ b_potential_dot_dot_acoustic)
+
+ use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
+ xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
+ station_name,network_name,adj_source_file,nrec_local,number_receiver_global, &
+ pm1_source_encoding,nsources_local,USE_FORCE_POINT_SOURCE, &
+ USE_RICKER_TIME_FUNCTION,factor_force_source
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSPEC_AB
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: kappastore
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! source
+ integer :: NSOURCES,myrank,it
+ integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
+ double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
+ double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,hdur_tiny,tshift_src
+ double precision :: dt,t0
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
+
+ double precision, external :: comp_source_time_function,comp_source_time_function_rickr,&
+ comp_source_time_function_gauss
+
+ logical, dimension(NSPEC_AB) :: ispec_is_acoustic
+
+!adjoint simulations
+ integer:: SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT
+ real(kind=CUSTOM_REAL),dimension(NGLOB_ADJOINT):: b_potential_dot_dot_acoustic
+
+! local parameters
+ double precision :: f0
+ double precision :: stf
+ real(kind=CUSTOM_REAL) stf_used,stf_used_total_all,time_source
+ integer :: isource,iglob,ispec,i,j,k
+
+! adjoint sources in SU format
+ integer,parameter :: nheader=240 ! 240 bytes
+
+! plotting source time function
+ if(PRINT_SOURCE_TIME_FUNCTION .and. .not. phase_is_inner ) then
+ ! initializes total
+ stf_used_total = 0.0_CUSTOM_REAL
+ endif
+
+! NOTE: adjoint sources and backward wavefield timing:
+! idea is to start with the backward field b_potential,.. at time (T)
+! and convolve with the adjoint field at time (T-t)
+!
+! backward/reconstructed wavefields:
+! time for b_potential( it ) would correspond to (NSTEP - it - 1 )*DT - t0
+! if we read in saved wavefields b_potential() before Newmark time scheme
+! (see sources for simulation_type 1 and seismograms)
+! since at the beginning of the time loop, the numerical Newmark time scheme updates
+! the wavefields, that is b_potential( it=1) would correspond to time (NSTEP -1 - 1)*DT - t0
+!
+! b_potential is now read in after Newmark time scheme:
+! we read the backward/reconstructed wavefield at the end of the first time loop,
+! such that b_potential(it=1) corresponds to -t0 + (NSTEP-1)*DT.
+! assuming that until that end the backward/reconstructed wavefield and adjoint fields
+! have a zero contribution to adjoint kernels.
+! thus the correct indexing is NSTEP - it + 1, instead of NSTEP - it
+!
+! adjoint wavefields:
+! since the adjoint source traces were derived from the seismograms,
+! it follows that for the adjoint wavefield, the time equivalent to ( T - t ) uses the time-reversed
+! adjoint source traces which start at -t0 and end at time (NSTEP-1)*DT - t0
+! for step it=1: (NSTEP -it + 1)*DT - t0 for backward wavefields corresponds to time T
+
! note: b_potential() is read in after Newmark time scheme, thus
! b_potential(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! thus indexing is NSTEP - it , instead of NSTEP - it - 1
@@ -470,7 +566,7 @@
if( myrank == 0 ) write(IOSTF,*) time_source,stf_used_total_all
endif
- end subroutine compute_add_sources_acoustic
+ end subroutine compute_add_sources_acoustic_bpwf
!
!=====================================================================
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_viscoelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_viscoelastic.f90 2013-06-08 12:12:16 UTC (rev 22193)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_viscoelastic.f90 2013-06-09 15:22:30 UTC (rev 22194)
@@ -23,16 +23,15 @@
! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
!
!=====================================================================
-
! for elastic solver
subroutine compute_add_sources_viscoelastic( NSPEC_AB,NGLOB_AB,accel, &
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
- ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ ispec_is_elastic,SIMULATION_TYPE,NSTEP, &
nrec,islice_selected_rec,ispec_selected_rec, &
- nadj_rec_local,adj_sourcearrays,b_accel, &
+ nadj_rec_local,adj_sourcearrays, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY)
use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
@@ -75,11 +74,10 @@
logical, dimension(NSPEC_AB) :: ispec_is_elastic
!adjoint simulations
- integer:: SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT
+ integer:: SIMULATION_TYPE,NSTEP
integer:: nrec
integer,dimension(nrec) :: islice_selected_rec,ispec_selected_rec
integer:: nadj_rec_local
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
logical :: ibool_read_adj_arrays
integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY
real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: &
@@ -349,10 +347,136 @@
endif ! nadj_rec_local
endif !adjoint
-! note: b_displ() is read in after Newmark time scheme, thus
-! b_displ(it=1) corresponds to -t0 + (NSTEP-1)*DT.
-! thus indexing is NSTEP - it , instead of NSTEP - it - 1
+ ! master prints out source time function to file
+ if(PRINT_SOURCE_TIME_FUNCTION .and. phase_is_inner) then
+ time_source = (it-1)*DT - t0
+ call sum_all_cr(stf_used_total,stf_used_total_all)
+ if( myrank == 0 ) write(IOSTF,*) time_source,stf_used_total_all
+ endif
+ ! for noise simulations
+ ! we have two loops indicated by phase_is_inner ("inner elements/points" or "boundary elements/points")
+ ! here, we only add those noise sources once, when we are calculating for boudanry points (phase_is_inner==.false.),
+ ! because boundary points are claculated first!
+ if( .not. phase_is_inner ) then
+ if ( NOISE_TOMOGRAPHY == 1 ) then
+ ! the first step of noise tomography is to use |S(\omega)|^2 as a point force source at one of the receivers.
+ ! hence, instead of a moment tensor 'sourcearrays', a 'noise_sourcearray' for a point force is needed.
+ ! furthermore, the CMTSOLUTION needs to be zero, i.e., no earthquakes.
+ ! now this must be manually set in DATA/CMTSOLUTION, by USERS.
+ call add_source_master_rec_noise(myrank,nrec, &
+ NSTEP,accel,noise_sourcearray, &
+ ibool,islice_selected_rec,ispec_selected_rec, &
+ it,irec_master_noise, &
+ NSPEC_AB,NGLOB_AB)
+ else if ( NOISE_TOMOGRAPHY == 2 ) then
+ ! second step of noise tomography, i.e., read the surface movie saved at every timestep
+ ! use the movie to drive the ensemble forward wavefield
+ call noise_read_add_surface_movie(NGLLSQUARE*num_free_surface_faces,accel, &
+ normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, &
+ ibool,noise_surface_movie,NSTEP-it+1,NSPEC_AB,NGLOB_AB, &
+ num_free_surface_faces,free_surface_ispec,free_surface_ijk, &
+ free_surface_jacobian2Dw)
+ ! be careful, since ensemble forward sources are reversals of generating wavefield "eta"
+ ! hence the "NSTEP-it+1", i.e., start reading from the last timestep
+ ! note the ensemble forward sources are generally distributed on the surface of the earth
+ ! that's to say, the ensemble forward source is kind of a surface force density, not a body force density
+ ! therefore, we must add it here, before applying the inverse of mass matrix
+ endif
+ endif
+
+
+ end subroutine compute_add_sources_viscoelastic
+!
+!=====================================================================
+! for elastic solver
+
+ subroutine compute_add_sources_viscoelastic_bpwf( NSPEC_AB,NGLOB_AB, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
+ hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
+ ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ b_accel,NOISE_TOMOGRAPHY)
+
+ use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
+ xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
+ station_name,network_name,adj_source_file, &
+ num_free_surface_faces,free_surface_ispec, &
+ free_surface_ijk,free_surface_jacobian2Dw, &
+ noise_sourcearray,irec_master_noise, &
+ normal_x_noise,normal_y_noise,normal_z_noise, &
+ mask_noise,noise_surface_movie, &
+ nrec_local,number_receiver_global, &
+ nsources_local,USE_FORCE_POINT_SOURCE, &
+ USE_RICKER_TIME_FUNCTION
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! source
+ integer :: NSOURCES,myrank,it
+ integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
+ double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,hdur_tiny,tshift_src
+ double precision :: dt,t0
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
+
+ double precision, external :: comp_source_time_function,comp_source_time_function_gauss,comp_source_time_function_rickr
+
+ logical, dimension(NSPEC_AB) :: ispec_is_elastic
+
+!adjoint simulations
+ integer:: SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
+ integer :: NOISE_TOMOGRAPHY
+
+! local parameters
+ double precision :: stf
+ real(kind=CUSTOM_REAL) stf_used,stf_used_total_all,time_source
+ integer :: isource,iglob,i,j,k,ispec
+
+! adjoint sources in SU format
+ integer,parameter :: nheader=240 ! 240 bytes
+
+! plotting source time function
+ if(PRINT_SOURCE_TIME_FUNCTION .and. .not. phase_is_inner ) then
+ ! initializes total
+ stf_used_total = 0.0_CUSTOM_REAL
+ endif
+
+! NOTE: adjoint sources and backward wavefield timing:
+! idea is to start with the backward field b_displ,.. at time (T)
+! and convolve with the adjoint field at time (T-t)
+!
+! backward/reconstructed wavefields:
+! time for b_displ( it ) would correspond to (NSTEP - it - 1 )*DT - t0
+! if we read in saved wavefields b_displ() before Newmark time scheme
+! (see sources for simulation_type 1 and seismograms)
+! since at the beginning of the time loop, the numerical Newmark time scheme updates
+! the wavefields, that is b_displ( it=1) would correspond to time (NSTEP -1 - 1)*DT - t0
+!
+! b_displ is now read in after Newmark time scheme:
+! we read the backward/reconstructed wavefield at the end of the first time loop,
+! such that b_displ(it=1) corresponds to -t0 + (NSTEP-1)*DT.
+! assuming that until that end the backward/reconstructed wavefield and adjoint fields
+! have a zero contribution to adjoint kernels.
+! thus the correct indexing is NSTEP - it + 1, instead of NSTEP - it
+!
+! adjoint wavefields:
+! since the adjoint source traces were derived from the seismograms,
+! it follows that for the adjoint wavefield, the time equivalent to ( T - t ) uses the time-reversed
+! adjoint source traces which start at -t0 and end at time (NSTEP-1)*DT - t0
+! for step it=1: (NSTEP -it + 1)*DT - t0 for backward wavefields corresponds to time T
+
! adjoint simulations
if (SIMULATION_TYPE == 3 .and. NOISE_TOMOGRAPHY == 0 .and. nsources_local > 0) then
@@ -450,30 +574,7 @@
! here, we only add those noise sources once, when we are calculating for boudanry points (phase_is_inner==.false.),
! because boundary points are claculated first!
if( .not. phase_is_inner ) then
- if ( NOISE_TOMOGRAPHY == 1 ) then
- ! the first step of noise tomography is to use |S(\omega)|^2 as a point force source at one of the receivers.
- ! hence, instead of a moment tensor 'sourcearrays', a 'noise_sourcearray' for a point force is needed.
- ! furthermore, the CMTSOLUTION needs to be zero, i.e., no earthquakes.
- ! now this must be manually set in DATA/CMTSOLUTION, by USERS.
- call add_source_master_rec_noise(myrank,nrec, &
- NSTEP,accel,noise_sourcearray, &
- ibool,islice_selected_rec,ispec_selected_rec, &
- it,irec_master_noise, &
- NSPEC_AB,NGLOB_AB)
- else if ( NOISE_TOMOGRAPHY == 2 ) then
- ! second step of noise tomography, i.e., read the surface movie saved at every timestep
- ! use the movie to drive the ensemble forward wavefield
- call noise_read_add_surface_movie(NGLLSQUARE*num_free_surface_faces,accel, &
- normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, &
- ibool,noise_surface_movie,NSTEP-it+1,NSPEC_AB,NGLOB_AB, &
- num_free_surface_faces,free_surface_ispec,free_surface_ijk, &
- free_surface_jacobian2Dw)
- ! be careful, since ensemble forward sources are reversals of generating wavefield "eta"
- ! hence the "NSTEP-it+1", i.e., start reading from the last timestep
- ! note the ensemble forward sources are generally distributed on the surface of the earth
- ! that's to say, the ensemble forward source is kind of a surface force density, not a body force density
- ! therefore, we must add it here, before applying the inverse of mass matrix
- else if ( NOISE_TOMOGRAPHY == 3 ) then
+ if ( NOISE_TOMOGRAPHY == 3 ) then
! third step of noise tomography, i.e., read the surface movie saved at every timestep
! use the movie to reconstruct the ensemble forward wavefield
! the ensemble adjoint wavefield is done as usual
@@ -487,7 +588,7 @@
endif
- end subroutine compute_add_sources_viscoelastic
+ end subroutine compute_add_sources_viscoelastic_bpwf
!
!=====================================================================
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90 2013-06-08 12:12:16 UTC (rev 22193)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90 2013-06-09 15:22:30 UTC (rev 22194)
@@ -142,8 +142,7 @@
ibool,rmassx,rmassy,rmassz, &
rmass_ocean_load,accel, &
free_surface_normal,free_surface_ijk,free_surface_ispec, &
- num_free_surface_faces,SIMULATION_TYPE, &
- NGLOB_ADJOINT,b_accel)
+ num_free_surface_faces)
! updates acceleration with ocean load term:
! approximates ocean-bottom continuity of pressure & displacement for longer period waves (> ~20s ),
@@ -167,10 +166,6 @@
integer :: free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces)
integer :: free_surface_ispec(num_free_surface_faces)
- ! adjoint simulations
- integer :: SIMULATION_TYPE,NGLOB_ADJOINT
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
-
! local parameters
real(kind=CUSTOM_REAL) :: nx,ny,nz
real(kind=CUSTOM_REAL) :: force_normal_comp
@@ -217,6 +212,85 @@
accel(3,iglob) = accel(3,iglob) &
+ (rmass_ocean_load(iglob) - rmassz(iglob)) * force_normal_comp * nz
+ ! done with this point
+ updated_dof_ocean_load(iglob) = .true.
+
+ endif
+
+ enddo ! igll
+ enddo ! iface
+
+ end subroutine compute_coupling_ocean
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine compute_coupling_ocean_bpwf(NSPEC_AB,NGLOB_AB, &
+ ibool,rmassx,rmassy,rmassz,rmass_ocean_load, &
+ free_surface_normal,free_surface_ijk,free_surface_ispec, &
+ num_free_surface_faces,SIMULATION_TYPE, &
+ NGLOB_ADJOINT,b_accel)
+
+! updates acceleration with ocean load term:
+! approximates ocean-bottom continuity of pressure & displacement for longer period waves (> ~20s ),
+! assuming incompressible fluid column above bathymetry ocean bottom
+
+ implicit none
+
+ include 'constants.h'
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+ real(kind=CUSTOM_REAL),dimension(NGLOB_AB),intent(in) :: rmassx,rmassy,rmassz
+ real(kind=CUSTOM_REAL),dimension(NGLOB_AB),intent(in) :: rmass_ocean_load
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: ibool
+
+ ! free surface
+ integer :: num_free_surface_faces
+ real(kind=CUSTOM_REAL) :: free_surface_normal(NDIM,NGLLSQUARE,num_free_surface_faces)
+ integer :: free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces)
+ integer :: free_surface_ispec(num_free_surface_faces)
+
+ ! adjoint simulations
+ integer :: SIMULATION_TYPE,NGLOB_ADJOINT
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
+
+! local parameters
+ real(kind=CUSTOM_REAL) :: nx,ny,nz
+ integer :: i,j,k,ispec,iglob
+ integer :: igll,iface
+ logical,dimension(NGLOB_AB) :: updated_dof_ocean_load
+ ! adjoint locals
+ real(kind=CUSTOM_REAL) :: b_force_normal_comp
+
+ ! initialize the updates
+ updated_dof_ocean_load(:) = .false.
+
+ ! for surface elements exactly at the top of the model (ocean bottom)
+ do iface = 1,num_free_surface_faces
+
+ ispec = free_surface_ispec(iface)
+ do igll = 1, NGLLSQUARE
+ i = free_surface_ijk(1,igll,iface)
+ j = free_surface_ijk(2,igll,iface)
+ k = free_surface_ijk(3,igll,iface)
+
+ ! get global point number
+ iglob = ibool(i,j,k,ispec)
+
+ ! only update once
+ if(.not. updated_dof_ocean_load(iglob)) then
+
+ ! get normal
+ nx = free_surface_normal(1,igll,iface)
+ ny = free_surface_normal(2,igll,iface)
+ nz = free_surface_normal(3,igll,iface)
+
+ ! make updated component of right-hand side
+ ! we divide by rmass() which is 1 / M
+ ! we use the total force which includes the Coriolis term above
+
! adjoint simulations
if (SIMULATION_TYPE == 3) then
b_force_normal_comp = b_accel(1,iglob)*nx / rmassx(iglob) &
@@ -239,5 +313,6 @@
enddo ! igll
enddo ! iface
- end subroutine compute_coupling_ocean
+ end subroutine compute_coupling_ocean_bpwf
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-06-08 12:12:16 UTC (rev 22193)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-06-09 15:22:30 UTC (rev 22194)
@@ -72,13 +72,6 @@
ibool,free_surface_ijk,free_surface_ispec, &
num_free_surface_faces,ispec_is_acoustic)
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) &
- call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_ADJOINT,STACEY_INSTEAD_OF_FREE_SURFACE, &
- b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
- ibool,free_surface_ijk,free_surface_ispec, &
- num_free_surface_faces,ispec_is_acoustic)
-
! distinguishes two runs: for elements on MPI interfaces, and elements within the partitions
do iphase=1,2
@@ -113,32 +106,6 @@
.false.,potential_dot_dot_acoustic_interface)
endif
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- if(USE_DEVILLE_PRODUCTS) then
- ! uses Deville (2002) optimizations
- call compute_forces_acoustic_Dev(iphase,NSPEC_ADJOINT,NGLOB_ADJOINT, &
- b_potential_acoustic,b_potential_dot_dot_acoustic, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D, &
- rhostore,jacobian,ibool, &
- num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
- phase_ispec_inner_acoustic)
- else
- call compute_forces_acoustic_noDev(iphase,NSPEC_ADJOINT,NGLOB_ADJOINT, &
- b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- hprime_xx,hprime_yy,hprime_zz, &
- hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- rhostore,jacobian,ibool,deltat, &
- num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
- phase_ispec_inner_acoustic,ELASTIC_SIMULATION,&
- .true.,potential_dot_dot_acoustic_interface)
- endif
- endif
-
! ! Stacey absorbing boundary conditions
if(STACEY_ABSORBING_CONDITIONS) then
call compute_stacey_acoustic(NSPEC_AB,NGLOB_AB, &
@@ -146,8 +113,7 @@
ibool,ispec_is_inner,phase_is_inner, &
abs_boundary_jacobian2Dw,abs_boundary_ijk,abs_boundary_ispec, &
num_abs_boundary_faces,rhostore,kappastore,ispec_is_acoustic, &
- SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it,NGLOB_ADJOINT, &
- b_potential_dot_dot_acoustic,b_reclen_potential, &
+ SIMULATION_TYPE,SAVE_FORWARD,it,b_reclen_potential, &
b_absorb_potential,b_num_abs_boundary_faces)
endif
@@ -170,7 +136,7 @@
! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
! adjoint definition: \partial_t^2 \bfs^\dagger=-\frac{1}{\rho}\bfnabla\phi^\dagger
call compute_coupling_acoustic_el(NSPEC_AB,NGLOB_AB, &
- ibool,-accel_adj_coupling,potential_dot_dot_acoustic, &
+ ibool,-accel,potential_dot_dot_acoustic, &
num_coupling_ac_el_faces, &
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
@@ -179,17 +145,6 @@
PML_CONDITIONS,spec_to_CPML,is_CPML,&
potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ)
endif
- ! adjoint/kernel simulations
- if( SIMULATION_TYPE == 3 ) &
- call compute_coupling_acoustic_el(NSPEC_ADJOINT,NGLOB_ADJOINT, &
- ibool,b_displ,b_potential_dot_dot_acoustic, &
- num_coupling_ac_el_faces, &
- coupling_ac_el_ispec,coupling_ac_el_ijk, &
- coupling_ac_el_normal, &
- coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,spec_to_CPML,is_CPML,&
- potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ)
endif
endif
@@ -208,8 +163,6 @@
else
stop 'not implemented yet'
endif
- if( SIMULATION_TYPE == 3 ) &
- stop 'not implemented yet'
endif
endif
@@ -220,10 +173,9 @@
xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
- SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ SIMULATION_TYPE,NSTEP, &
nrec,islice_selected_rec,ispec_selected_rec, &
- nadj_rec_local,adj_sourcearrays,b_potential_dot_dot_acoustic, &
- NTSTEP_BETWEEN_READ_ADJSRC)
+ nadj_rec_local,adj_sourcearrays,NTSTEP_BETWEEN_READ_ADJSRC)
! assemble all the contributions between slices using MPI
if( phase_is_inner .eqv. .false. ) then
@@ -235,16 +187,6 @@
my_neighbours_ext_mesh, &
request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh)
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- call assemble_MPI_scalar_ext_mesh_s(NPROC,NGLOB_ADJOINT,b_potential_dot_dot_acoustic, &
- b_buffer_send_scalar_ext_mesh,b_buffer_recv_scalar_ext_mesh, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
- my_neighbours_ext_mesh, &
- b_request_send_scalar_ext_mesh,b_request_recv_scalar_ext_mesh)
- endif
-
else
! waits for send/receive requests to be completed and assembles values
@@ -254,14 +196,6 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh)
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- call assemble_MPI_scalar_ext_mesh_w(NPROC,NGLOB_ADJOINT,b_potential_dot_dot_acoustic, &
- b_buffer_recv_scalar_ext_mesh,num_interfaces_ext_mesh,&
- max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
- b_request_send_scalar_ext_mesh,b_request_recv_scalar_ext_mesh)
- endif
endif !phase_is_inner
enddo
@@ -276,10 +210,6 @@
endif
endif
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) &
- b_potential_dot_dot_acoustic(:) = b_potential_dot_dot_acoustic(:) * rmass_acoustic(:)
-
! The outer boundary condition to use for PML elements in fluid layers is Neumann for the potential
! because we need Dirichlet conditions for the displacement vector, which means Neumann for the potential.
! Thus, there is nothing to enforce explicitly here.
@@ -336,10 +266,6 @@
! corrector
potential_dot_acoustic(:) = potential_dot_acoustic(:) + deltatover2*potential_dot_dot_acoustic(:)
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) &
- b_potential_dot_acoustic(:) = b_potential_dot_acoustic(:) + b_deltatover2*b_potential_dot_dot_acoustic(:)
-
! enforces free surface (zeroes potentials at free surface)
call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_AB,STACEY_INSTEAD_OF_FREE_SURFACE, &
potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
@@ -352,13 +278,6 @@
+ deltatsqover2 * potential_dot_dot_acoustic(:)
endif
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) &
- call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_ADJOINT,STACEY_INSTEAD_OF_FREE_SURFACE, &
- b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
- ibool,free_surface_ijk,free_surface_ispec, &
- num_free_surface_faces,ispec_is_acoustic)
-
if(PML_CONDITIONS)then
if(SIMULATION_TYPE == 1 .and. SAVE_FORWARD)then
if(nglob_interface_PML_acoustic > 0)then
@@ -371,63 +290,200 @@
end subroutine compute_forces_acoustic
!
-!-------------------------------------------------------------------------------------------------
+!=====================================================================
+
+! acoustic solver
+
+! in case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette,
+! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential
+! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002).
!
+! This permits acoustic-elastic coupling based on a non-iterative time scheme.
+! Displacement is then:
+! u = grad(Chi) / rho
+! Velocity is then:
+! v = grad(Chi_dot) / rho
+! (Chi_dot being the time derivative of Chi)
+! and pressure is:
+! p = - Chi_dot_dot
+! (Chi_dot_dot being the time second derivative of Chi).
+!
+! The source in an acoustic element is a pressure source.
+!
+! First-order acoustic-acoustic discontinuities are also handled automatically
+! because pressure is continuous at such an interface, therefore Chi_dot_dot
+! is continuous, therefore Chi is also continuous, which is consistent with
+! the spectral-element basis functions and with the assembling process.
+! This is the reason why a simple displacement potential u = grad(Chi) would
+! not work because it would be discontinuous at such an interface and would
+! therefore not be consistent with the basis functions.
-subroutine acoustic_enforce_free_surface(NSPEC_AB,NGLOB_AB,STACEY_INSTEAD_OF_FREE_SURFACE, &
- potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
- ibool,free_surface_ijk,free_surface_ispec, &
- num_free_surface_faces,ispec_is_acoustic)
+
+subroutine compute_forces_acoustic_bpwf()
+
+ use specfem_par
+ use specfem_par_acoustic
+ use specfem_par_elastic
+ use specfem_par_poroelastic
+ use pml_par,only: spec_to_CPML,is_CPML,rmemory_dpotential_dxl,rmemory_dpotential_dyl,rmemory_dpotential_dzl,&
+ rmemory_potential_acoustic,rmemory_coupling_ac_el_displ,nglob_interface_PML_acoustic,&
+ b_PML_potential,b_reclen_PML_potential
implicit none
- include 'constants.h'
- integer :: NSPEC_AB,NGLOB_AB
- logical :: STACEY_INSTEAD_OF_FREE_SURFACE
+ ! local parameters
+ integer:: iphase
+ logical:: phase_is_inner
-! acoustic potentials
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: &
- potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) &
+ call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_ADJOINT,STACEY_INSTEAD_OF_FREE_SURFACE, &
+ b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
+ ibool,free_surface_ijk,free_surface_ispec, &
+ num_free_surface_faces,ispec_is_acoustic)
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ ! distinguishes two runs: for elements on MPI interfaces, and elements within the partitions
+ do iphase=1,2
-! free surface
- integer :: num_free_surface_faces
- integer :: free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces)
- integer :: free_surface_ispec(num_free_surface_faces)
+ !first for points on MPI interfaces, thus outer elements
+ if( iphase == 1 ) then
+ phase_is_inner = .false.
+ else
+ phase_is_inner = .true.
+ endif
- logical, dimension(NSPEC_AB) :: ispec_is_acoustic
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) then
+ if(USE_DEVILLE_PRODUCTS) then
+ ! uses Deville (2002) optimizations
+ call compute_forces_acoustic_Dev(iphase,NSPEC_ADJOINT,NGLOB_ADJOINT, &
+ b_potential_acoustic,b_potential_dot_dot_acoustic, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D, &
+ rhostore,jacobian,ibool, &
+ num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
+ phase_ispec_inner_acoustic)
+ else
+ call compute_forces_acoustic_noDev(iphase,NSPEC_ADJOINT,NGLOB_ADJOINT, &
+ b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ rhostore,jacobian,ibool,deltat, &
+ num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
+ phase_ispec_inner_acoustic,ELASTIC_SIMULATION,&
+ .true.,potential_dot_dot_acoustic_interface)
+ endif
+ endif
-! local parameters
- integer :: iface,igll,i,j,k,ispec,iglob
+ ! ! Stacey absorbing boundary conditions
+ if(STACEY_ABSORBING_CONDITIONS) then
+ call compute_stacey_acoustic_bpwf(NSPEC_AB, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ abs_boundary_ijk,abs_boundary_ispec, &
+ num_abs_boundary_faces,ispec_is_acoustic,&
+ SIMULATION_TYPE,NSTEP,it,NGLOB_ADJOINT, &
+ b_potential_dot_dot_acoustic,b_reclen_potential, &
+ b_absorb_potential,b_num_abs_boundary_faces)
+ endif
- ! checks if free surface became an absorbing boundary
- if( STACEY_INSTEAD_OF_FREE_SURFACE ) return
+ ! elastic coupling
+ if(ELASTIC_SIMULATION ) then
+ if( num_coupling_ac_el_faces > 0 ) then
+ ! adjoint/kernel simulations
+ if( SIMULATION_TYPE == 3 ) &
+ call compute_coupling_acoustic_el(NSPEC_ADJOINT,NGLOB_ADJOINT, &
+ ibool,b_displ,b_potential_dot_dot_acoustic, &
+ num_coupling_ac_el_faces, &
+ coupling_ac_el_ispec,coupling_ac_el_ijk, &
+ coupling_ac_el_normal, &
+ coupling_ac_el_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner,&
+ PML_CONDITIONS,spec_to_CPML,is_CPML,&
+ potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ)
+ endif
+ endif
-! enforce potentials to be zero at surface
- do iface = 1, num_free_surface_faces
+! poroelastic coupling
+ if(POROELASTIC_SIMULATION ) then
+ if( num_coupling_ac_po_faces > 0 ) then
+ if( SIMULATION_TYPE == 3 ) &
+ stop 'not implemented yet'
+ endif
+ endif
- ispec = free_surface_ispec(iface)
+ ! sources
+ call compute_add_sources_acoustic_bpwf(NSPEC_AB, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
+ xi_source,eta_source,gamma_source, &
+ hdur,hdur_gaussian,tshift_src,dt,t0, &
+ sourcearrays,kappastore,ispec_is_acoustic,&
+ SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ b_potential_dot_dot_acoustic)
- if( ispec_is_acoustic(ispec) ) then
+ ! assemble all the contributions between slices using MPI
+ if( phase_is_inner .eqv. .false. ) then
+ ! sends b_potential_dot_dot_acoustic values to corresponding MPI interface neighbors (non-blocking)
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) then
+ call assemble_MPI_scalar_ext_mesh_s(NPROC,NGLOB_ADJOINT,b_potential_dot_dot_acoustic, &
+ b_buffer_send_scalar_ext_mesh,b_buffer_recv_scalar_ext_mesh, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
+ my_neighbours_ext_mesh, &
+ b_request_send_scalar_ext_mesh,b_request_recv_scalar_ext_mesh)
+ endif
- do igll = 1, NGLLSQUARE
- i = free_surface_ijk(1,igll,iface)
- j = free_surface_ijk(2,igll,iface)
- k = free_surface_ijk(3,igll,iface)
- iglob = ibool(i,j,k,ispec)
+ else
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) then
+ call assemble_MPI_scalar_ext_mesh_w(NPROC,NGLOB_ADJOINT,b_potential_dot_dot_acoustic, &
+ b_buffer_recv_scalar_ext_mesh,num_interfaces_ext_mesh,&
+ max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ b_request_send_scalar_ext_mesh,b_request_recv_scalar_ext_mesh)
+ endif
+ endif !phase_is_inner
- ! sets potentials to zero
- potential_acoustic(iglob) = 0._CUSTOM_REAL
- potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
- potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
- enddo
- endif
-
enddo
-end subroutine acoustic_enforce_free_surface
+ ! divides pressure with mass matrix
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) &
+ b_potential_dot_dot_acoustic(:) = b_potential_dot_dot_acoustic(:) * rmass_acoustic(:)
+! update velocity
+! note: Newmark finite-difference time scheme with acoustic domains:
+! (see e.g. Hughes, 1987; Chaljub et al., 2003)
!
+! chi(t+delta_t) = chi(t) + delta_t chi_dot(t) + 1/2 delta_t**2 chi_dot_dot(t)
+! chi_dot(t+delta_t) = chi_dot(t) + 1/2 delta_t chi_dot_dot(t) + 1/2 DELTA_T CHI_DOT_DOT( T + DELTA_T )
+! chi_dot_dot(t+delta_t) = 1/M_acoustic( -K_acoustic chi(t+delta) + B_acoustic u(t+delta_t) + f(t+delta_t) )
+!
+! where
+! chi, chi_dot, chi_dot_dot are acoustic (fluid) potentials ( dotted with respect to time)
+! M is mass matrix, K stiffness matrix and B boundary term
+! f denotes a source term
+!
+! corrector:
+! updates the chi_dot term which requires chi_dot_dot(t+delta)
+ ! corrector
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) &
+ b_potential_dot_acoustic(:) = b_potential_dot_acoustic(:) + b_deltatover2*b_potential_dot_dot_acoustic(:)
+
+! enforces free surface (zeroes potentials at free surface)
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) &
+ call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_ADJOINT,STACEY_INSTEAD_OF_FREE_SURFACE, &
+ b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
+ ibool,free_surface_ijk,free_surface_ispec, &
+ num_free_surface_faces,ispec_is_acoustic)
+
+end subroutine compute_forces_acoustic_bpwf
+!
!-------------------------------------------------------------------------------------------------
!
subroutine compute_forces_acoustic_GPU()
@@ -464,7 +520,7 @@
! ! Stacey absorbing boundary conditions
if(STACEY_ABSORBING_CONDITIONS) then
call compute_stacey_acoustic_GPU(phase_is_inner,num_abs_boundary_faces,&
- SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it,NGLOB_ADJOINT, &
+ SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it, &
b_reclen_potential,b_absorb_potential, &
b_num_abs_boundary_faces,Mesh_pointer)
endif
@@ -594,4 +650,60 @@
call acoustic_enforce_free_surf_cuda(Mesh_pointer,STACEY_INSTEAD_OF_FREE_SURFACE)
end subroutine compute_forces_acoustic_GPU
+!
+!-------------------------------------------------------------------------------------------------
+!
+subroutine acoustic_enforce_free_surface(NSPEC_AB,NGLOB_AB,STACEY_INSTEAD_OF_FREE_SURFACE, &
+ potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
+ ibool,free_surface_ijk,free_surface_ispec, &
+ num_free_surface_faces,ispec_is_acoustic)
+ implicit none
+ include 'constants.h'
+
+ integer :: NSPEC_AB,NGLOB_AB
+ logical :: STACEY_INSTEAD_OF_FREE_SURFACE
+
+! acoustic potentials
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: &
+ potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+! free surface
+ integer :: num_free_surface_faces
+ integer :: free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces)
+ integer :: free_surface_ispec(num_free_surface_faces)
+
+ logical, dimension(NSPEC_AB) :: ispec_is_acoustic
+
+! local parameters
+ integer :: iface,igll,i,j,k,ispec,iglob
+
+ ! checks if free surface became an absorbing boundary
+ if( STACEY_INSTEAD_OF_FREE_SURFACE ) return
+
+! enforce potentials to be zero at surface
+ do iface = 1, num_free_surface_faces
+
+ ispec = free_surface_ispec(iface)
+
+ if( ispec_is_acoustic(ispec) ) then
+
+ do igll = 1, NGLLSQUARE
+ i = free_surface_ijk(1,igll,iface)
+ j = free_surface_ijk(2,igll,iface)
+ k = free_surface_ijk(3,igll,iface)
+ iglob = ibool(i,j,k,ispec)
+
+ ! sets potentials to zero
+ potential_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ enddo
+ endif
+
+ enddo
+
+end subroutine acoustic_enforce_free_surface
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-06-08 12:12:16 UTC (rev 22193)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-06-09 15:22:30 UTC (rev 22194)
@@ -52,15 +52,11 @@
phase_is_inner = .true.
endif
-
-! elastic term
+ ! elastic term
if(USE_DEVILLE_PRODUCTS) then
! uses Deville (2002) optimizations
call compute_forces_viscoelastic_Dev_sim1(iphase)
- ! adjoint simulations: backward/reconstructed wavefield
- if( SIMULATION_TYPE == 3 ) call compute_forces_viscoelastic_Dev_sim3(iphase)
-
else
! no optimizations used
call compute_forces_viscoelastic_noDev(iphase,NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
@@ -90,36 +86,6 @@
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
phase_ispec_inner_elastic,.false.)
- ! adjoint simulations: backward/reconstructed wavefield
- if( SIMULATION_TYPE == 3 ) &
- call compute_forces_viscoelastic_noDev(iphase,NSPEC_AB,NGLOB_AB, &
- b_displ,b_veloc,b_accel, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- hprime_xx,hprime_yy,hprime_zz, &
- hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- kappastore,mustore,jacobian,ibool, &
- ATTENUATION,deltat,PML_CONDITIONS, &
- one_minus_sum_beta,factor_common, &
- one_minus_sum_beta_kappa,factor_common_kappa, &
- b_alphaval,b_betaval,b_gammaval, &
- NSPEC_ATTENUATION_AB,NSPEC_ATTENUATION_AB_kappa, &
- b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
- b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
- b_epsilondev_xz,b_epsilondev_yz,b_epsilon_trace_over_3, &
- ANISOTROPY,NSPEC_ANISO, &
- c11store,c12store,c13store,c14store,c15store,c16store, &
- c22store,c23store,c24store,c25store,c26store,c33store, &
- c34store,c35store,c36store,c44store,c45store,c46store, &
- c55store,c56store,c66store, &
- SIMULATION_TYPE,COMPUTE_AND_STORE_STRAIN,NSPEC_STRAIN_ONLY, &
- NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ADJOINT, &
- is_moho_top,is_moho_bot, &
- b_dsdx_top,b_dsdx_bot, &
- ispec2D_moho_top,ispec2D_moho_bot, &
- num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,.true.)
-
endif
@@ -132,7 +98,7 @@
num_abs_boundary_faces, &
veloc,rho_vp,rho_vs, &
ispec_is_elastic,SIMULATION_TYPE,SAVE_FORWARD, &
- NSTEP,it,NGLOB_ADJOINT,b_accel, &
+ NSTEP,it, &
b_num_abs_boundary_faces,b_reclen_field,b_absorb_field,&
it_dsm,Veloc_dsm_boundary,Tract_dsm_boundary)
endif
@@ -155,7 +121,7 @@
! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
! adoint definition: pressure^\dagger=potential^\dagger
call compute_coupling_viscoelastic_ac(NSPEC_AB,NGLOB_AB, &
- ibool,accel,-potential_acoustic_adj_coupling, &
+ ibool,accel,-potential_acoustic, &
num_coupling_ac_el_faces, &
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
@@ -164,17 +130,6 @@
PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface)
endif
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) &
- call compute_coupling_viscoelastic_ac(NSPEC_ADJOINT,NGLOB_ADJOINT, &
- ibool,b_accel,b_potential_dot_dot_acoustic, &
- num_coupling_ac_el_faces, &
- coupling_ac_el_ispec,coupling_ac_el_ijk, &
- coupling_ac_el_normal, &
- coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface)
-
endif ! num_coupling_ac_el_faces
endif
@@ -207,9 +162,9 @@
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
- ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ ispec_is_elastic,SIMULATION_TYPE,NSTEP, &
nrec,islice_selected_rec,ispec_selected_rec, &
- nadj_rec_local,adj_sourcearrays,b_accel, &
+ nadj_rec_local,adj_sourcearrays, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY)
! assemble all the contributions between slices using MPI
@@ -222,16 +177,6 @@
my_neighbours_ext_mesh, &
request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- call assemble_MPI_vector_ext_mesh_s(NPROC,NGLOB_ADJOINT,b_accel, &
- b_buffer_send_vector_ext_mesh,b_buffer_recv_vector_ext_mesh, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
- my_neighbours_ext_mesh, &
- b_request_send_vector_ext_mesh,b_request_recv_vector_ext_mesh)
- endif !adjoint
-
else
! waits for send/receive requests to be completed and assembles values
call assemble_MPI_vector_ext_mesh_w_ordered(NPROC,NGLOB_AB,accel, &
@@ -240,16 +185,6 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
request_send_vector_ext_mesh,request_recv_vector_ext_mesh, &
my_neighbours_ext_mesh,myrank)
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- call assemble_MPI_vector_ext_mesh_w_ordered(NPROC,NGLOB_ADJOINT,b_accel, &
- b_buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh,&
- max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
- b_request_send_vector_ext_mesh,b_request_recv_vector_ext_mesh, &
- my_neighbours_ext_mesh,myrank)
- endif !adjoint
-
endif
enddo
@@ -264,12 +199,6 @@
accel(1,:) = accel(1,:)*rmassx(:)
accel(2,:) = accel(2,:)*rmassy(:)
accel(3,:) = accel(3,:)*rmassz(:)
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_accel(1,:) = b_accel(1,:)*rmassx(:)
- b_accel(2,:) = b_accel(2,:)*rmassy(:)
- b_accel(3,:) = b_accel(3,:)*rmassz(:)
- endif !adjoint
! updates acceleration with ocean load term
if(APPROXIMATE_OCEAN_LOAD) then
@@ -277,8 +206,7 @@
ibool,rmassx,rmassy,rmassz, &
rmass_ocean_load,accel, &
free_surface_normal,free_surface_ijk,free_surface_ispec, &
- num_free_surface_faces,SIMULATION_TYPE, &
- NGLOB_ADJOINT,b_accel)
+ num_free_surface_faces)
endif
! C-PML boundary
@@ -325,8 +253,6 @@
! corrector:
! updates the velocity term which requires a(t+delta)
veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) b_veloc(:,:) = b_veloc(:,:) + b_deltatover2*b_accel(:,:)
if(PML_CONDITIONS)then
if(SIMULATION_TYPE == 1 .and. SAVE_FORWARD)then
@@ -339,6 +265,187 @@
end subroutine compute_forces_viscoelastic
!
+!=====================================================================
+
+! elastic solver for back
+
+subroutine compute_forces_viscoelastic_bpwf()
+
+ use specfem_par
+ use specfem_par_acoustic
+ use specfem_par_elastic
+ use specfem_par_poroelastic
+ use pml_par
+ use fault_solver_dynamic, only : bc_dynflt_set3d_all,SIMULATION_TYPE_DYN
+ use fault_solver_kinematic, only : bc_kinflt_set_all,SIMULATION_TYPE_KIN
+
+ implicit none
+
+ integer:: iphase
+ logical:: phase_is_inner
+
+! distinguishes two runs: for points on MPI interfaces, and points within the partitions
+ do iphase=1,2
+
+ !first for points on MPI interfaces
+ if( iphase == 1 ) then
+ phase_is_inner = .false.
+ else
+ phase_is_inner = .true.
+ endif
+
+
+! elastic term
+ if(USE_DEVILLE_PRODUCTS) then
+
+ ! adjoint simulations: backward/reconstructed wavefield
+ if( SIMULATION_TYPE == 3 ) call compute_forces_viscoelastic_Dev_sim3(iphase)
+
+ else
+ ! no optimizations used
+ ! adjoint simulations: backward/reconstructed wavefield
+ if( SIMULATION_TYPE == 3 ) &
+ call compute_forces_viscoelastic_noDev(iphase,NSPEC_AB,NGLOB_AB, &
+ b_displ,b_veloc,b_accel, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ kappastore,mustore,jacobian,ibool, &
+ ATTENUATION,deltat,PML_CONDITIONS, &
+ one_minus_sum_beta,factor_common, &
+ one_minus_sum_beta_kappa,factor_common_kappa, &
+ b_alphaval,b_betaval,b_gammaval, &
+ NSPEC_ATTENUATION_AB,NSPEC_ATTENUATION_AB_kappa, &
+ b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
+ b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
+ b_epsilondev_xz,b_epsilondev_yz,b_epsilon_trace_over_3, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store, &
+ SIMULATION_TYPE,COMPUTE_AND_STORE_STRAIN,NSPEC_STRAIN_ONLY, &
+ NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ADJOINT, &
+ is_moho_top,is_moho_bot, &
+ b_dsdx_top,b_dsdx_bot, &
+ ispec2D_moho_top,ispec2D_moho_bot, &
+ num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
+ phase_ispec_inner_elastic,.true.)
+
+ endif
+
+
+! adds elastic absorbing boundary term to acceleration (Stacey conditions)
+ if( STACEY_ABSORBING_CONDITIONS ) then
+ call compute_stacey_viscoelastic_bpwf(NSPEC_AB, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ abs_boundary_ijk,abs_boundary_ispec, &
+ num_abs_boundary_faces, &
+ ispec_is_elastic,SIMULATION_TYPE, &
+ NSTEP,it,NGLOB_ADJOINT,b_accel, &
+ b_num_abs_boundary_faces,b_reclen_field,b_absorb_field)
+ endif
+
+
+! acoustic coupling
+ if( ACOUSTIC_SIMULATION ) then
+ if( num_coupling_ac_el_faces > 0 ) then
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) &
+ call compute_coupling_viscoelastic_ac(NSPEC_ADJOINT,NGLOB_ADJOINT, &
+ ibool,b_accel,b_potential_dot_dot_acoustic, &
+ num_coupling_ac_el_faces, &
+ coupling_ac_el_ispec,coupling_ac_el_ijk, &
+ coupling_ac_el_normal, &
+ coupling_ac_el_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner,&
+ PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface)
+ endif ! num_coupling_ac_el_faces
+ endif
+
+
+! poroelastic coupling
+ if( POROELASTIC_SIMULATION ) then
+ stop 'poroelastic-elastic coupling error'
+ endif
+
+! adds source term (single-force/moment-tensor solution)
+ if (.not. OLD_TEST_TO_FIX_ONE_DAY) call compute_add_sources_viscoelastic_bpwf( NSPEC_AB,NGLOB_AB, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
+ hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
+ ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ b_accel,NOISE_TOMOGRAPHY)
+
+ ! assemble all the contributions between slices using MPI
+ if( phase_is_inner .eqv. .false. ) then
+ ! sends accel values to corresponding MPI interface neighbors
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) then
+ call assemble_MPI_vector_ext_mesh_s(NPROC,NGLOB_ADJOINT,b_accel, &
+ b_buffer_send_vector_ext_mesh,b_buffer_recv_vector_ext_mesh, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
+ my_neighbours_ext_mesh, &
+ b_request_send_vector_ext_mesh,b_request_recv_vector_ext_mesh)
+ endif !adjoint
+
+ else
+ ! waits for send/receive requests to be completed and assembles values
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) then
+ call assemble_MPI_vector_ext_mesh_w_ordered(NPROC,NGLOB_ADJOINT,b_accel, &
+ b_buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh,&
+ max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ b_request_send_vector_ext_mesh,b_request_recv_vector_ext_mesh, &
+ my_neighbours_ext_mesh,myrank)
+ endif !adjoint
+
+ endif
+
+ enddo
+
+ ! multiplies with inverse of mass matrix (note: rmass has been inverted already)
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ b_accel(1,:) = b_accel(1,:)*rmassx(:)
+ b_accel(2,:) = b_accel(2,:)*rmassy(:)
+ b_accel(3,:) = b_accel(3,:)*rmassz(:)
+ endif !adjoint
+
+! updates acceleration with ocean load term
+ if(APPROXIMATE_OCEAN_LOAD) then
+ call compute_coupling_ocean_bpwf(NSPEC_AB,NGLOB_AB, &
+ ibool,rmassx,rmassy,rmassz,rmass_ocean_load, &
+ free_surface_normal,free_surface_ijk,free_surface_ispec, &
+ num_free_surface_faces,SIMULATION_TYPE, &
+ NGLOB_ADJOINT,b_accel)
+ endif
+
+
+! updates velocities
+! Newmark finite-difference time scheme with elastic domains:
+! (see e.g. Hughes, 1987; Chaljub et al., 2003)
+!
+! u(t+delta_t) = u(t) + delta_t v(t) + 1/2 delta_t**2 a(t)
+! v(t+delta_t) = v(t) + 1/2 delta_t a(t) + 1/2 delta_t a(t+delta_t)
+! a(t+delta_t) = 1/M_elastic ( -K_elastic u(t+delta) + B_elastic chi_dot_dot(t+delta_t) + f( t+delta_t) )
+!
+! where
+! u, v, a are displacement,velocity & acceleration
+! M is mass matrix, K stiffness matrix and B boundary term for acoustic/elastic domains
+! f denotes a source term (acoustic/elastic)
+! chi_dot_dot is acoustic (fluid) potential ( dotted twice with respect to time)
+!
+! corrector:
+! updates the velocity term which requires a(t+delta)
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) b_veloc(:,:) = b_veloc(:,:) + b_deltatover2*b_accel(:,:)
+
+end subroutine compute_forces_viscoelastic_bpwf
+!
!-------------------------------------------------------------------------------------------------
!
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90 2013-06-08 12:12:16 UTC (rev 22193)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90 2013-06-09 15:22:30 UTC (rev 22194)
@@ -31,8 +31,7 @@
ibool,ispec_is_inner,phase_is_inner, &
abs_boundary_jacobian2Dw,abs_boundary_ijk,abs_boundary_ispec, &
num_abs_boundary_faces,rhostore,kappastore,ispec_is_acoustic,&
- SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it,NGLOB_ADJOINT, &
- b_potential_dot_dot_acoustic,b_reclen_potential, &
+ SIMULATION_TYPE,SAVE_FORWARD,it,b_reclen_potential, &
b_absorb_potential,b_num_abs_boundary_faces)
implicit none
@@ -60,10 +59,8 @@
integer :: abs_boundary_ispec(num_abs_boundary_faces)
! adjoint simulations
- integer:: SIMULATION_TYPE
- integer:: NSTEP,it,NGLOB_ADJOINT
+ integer:: SIMULATION_TYPE,it
integer:: b_num_abs_boundary_faces,b_reclen_potential
- real(kind=CUSTOM_REAL),dimension(NGLOB_ADJOINT) :: b_potential_dot_dot_acoustic
real(kind=CUSTOM_REAL),dimension(NGLLSQUARE,b_num_abs_boundary_faces):: b_absorb_potential
logical:: SAVE_FORWARD
@@ -75,20 +72,6 @@
! checks if anything to do
if( num_abs_boundary_faces == 0 ) return
- ! adjoint simulations:
- if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
- ! reads in absorbing boundary array when first phase is running
- if( phase_is_inner .eqv. .false. ) then
- ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newmark scheme
- ! uses fortran routine
- !read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
- !if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
- ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_potential')
- ! uses c routine for faster reading
- call read_abs(1,b_absorb_potential,b_reclen_potential,NSTEP-it+1)
- endif
- endif !adjoint
-
! absorbs absorbing-boundary surface using Sommerfeld condition (vanishing field in the outer-space)
do iface=1,num_abs_boundary_faces
@@ -121,12 +104,8 @@
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - absorbl
! adjoint simulations
- if(SIMULATION_TYPE == 3) then
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
- - b_absorb_potential(igll,iface)
- else if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD) then
+ if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD) then
b_absorb_potential(igll,iface) = absorbl
-
endif !adjoint
enddo
@@ -148,10 +127,102 @@
end subroutine compute_stacey_acoustic
!
!=====================================================================
+! for acoustic solver for back propagation wave field
+
+ subroutine compute_stacey_acoustic_bpwf(NSPEC_AB, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ abs_boundary_ijk,abs_boundary_ispec, &
+ num_abs_boundary_faces,ispec_is_acoustic,&
+ SIMULATION_TYPE,NSTEP,it,NGLOB_ADJOINT, &
+ b_potential_dot_dot_acoustic,b_reclen_potential, &
+ b_absorb_potential,b_num_abs_boundary_faces)
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSPEC_AB
+
+! potentials
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+ logical, dimension(NSPEC_AB) :: ispec_is_acoustic
+
+! absorbing boundary surface
+ integer :: num_abs_boundary_faces
+ integer :: abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces)
+ integer :: abs_boundary_ispec(num_abs_boundary_faces)
+
+! adjoint simulations
+ integer:: SIMULATION_TYPE
+ integer:: NSTEP,it,NGLOB_ADJOINT
+ integer:: b_num_abs_boundary_faces,b_reclen_potential
+ real(kind=CUSTOM_REAL),dimension(NGLOB_ADJOINT) :: b_potential_dot_dot_acoustic
+ real(kind=CUSTOM_REAL),dimension(NGLLSQUARE,b_num_abs_boundary_faces):: b_absorb_potential
+
+! local parameters
+ integer :: ispec,iglob,i,j,k,iface,igll
+ !integer:: reclen1,reclen2
+
+ ! checks if anything to do
+ if( num_abs_boundary_faces == 0 ) return
+
+ ! adjoint simulations:
+ if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
+ ! reads in absorbing boundary array when first phase is running
+ if( phase_is_inner .eqv. .false. ) then
+ ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newmark scheme
+ ! uses fortran routine
+ !read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
+ !if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
+ ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_potential')
+ ! uses c routine for faster reading
+ call read_abs(1,b_absorb_potential,b_reclen_potential,NSTEP-it+1)
+ endif
+ endif !adjoint
+
+ ! absorbs absorbing-boundary surface using Sommerfeld condition (vanishing field in the outer-space)
+ do iface=1,num_abs_boundary_faces
+
+ ispec = abs_boundary_ispec(iface)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ if( ispec_is_acoustic(ispec) ) then
+
+ ! reference gll points on boundary face
+ do igll = 1,NGLLSQUARE
+
+ ! gets local indices for GLL point
+ i = abs_boundary_ijk(1,igll,iface)
+ j = abs_boundary_ijk(2,igll,iface)
+ k = abs_boundary_ijk(3,igll,iface)
+
+ ! gets global index
+ iglob=ibool(i,j,k,ispec)
+
+ ! adjoint simulations
+ if(SIMULATION_TYPE == 3) then
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
+ - b_absorb_potential(igll,iface)
+ endif !adjoint
+
+ enddo
+ endif ! ispec_is_acoustic
+ endif ! ispec_is_inner
+ enddo ! num_abs_boundary_faces
+
+ end subroutine compute_stacey_acoustic_bpwf
+!
+!=====================================================================
! for acoustic solver on GPU
subroutine compute_stacey_acoustic_GPU(phase_is_inner,num_abs_boundary_faces,&
- SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it,NGLOB_ADJOINT, &
+ SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it, &
b_reclen_potential,b_absorb_potential, &
b_num_abs_boundary_faces,Mesh_pointer)
@@ -169,7 +240,7 @@
! adjoint simulations
integer:: SIMULATION_TYPE
- integer:: NSTEP,it,NGLOB_ADJOINT
+ integer:: NSTEP,it
integer:: b_num_abs_boundary_faces,b_reclen_potential
real(kind=CUSTOM_REAL),dimension(NGLLSQUARE,b_num_abs_boundary_faces):: b_absorb_potential
logical:: SAVE_FORWARD
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_viscoelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_viscoelastic.f90 2013-06-08 12:12:16 UTC (rev 22193)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_viscoelastic.f90 2013-06-09 15:22:30 UTC (rev 22194)
@@ -35,7 +35,7 @@
num_abs_boundary_faces, &
veloc,rho_vp,rho_vs, &
ispec_is_elastic,SIMULATION_TYPE,SAVE_FORWARD, &
- NSTEP,it,NGLOB_ADJOINT,b_accel, &
+ NSTEP,it, &
b_num_abs_boundary_faces,b_reclen_field,b_absorb_field, &
it_dsm,Veloc_dsm_boundary,Tract_dsm_boundary)
@@ -72,7 +72,6 @@
integer:: b_num_abs_boundary_faces,b_reclen_field
real(kind=CUSTOM_REAL),dimension(NDIM,NGLLSQUARE,b_num_abs_boundary_faces):: b_absorb_field
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
logical:: SAVE_FORWARD
! local parameters
@@ -171,9 +170,7 @@
accel(3,iglob) = accel(3,iglob) - tz*jacobianw
! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_accel(:,iglob) = b_accel(:,iglob) - b_absorb_field(:,igll,iface)
- else if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD) then
+ if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD) then
b_absorb_field(1,igll,iface) = tx*jacobianw
b_absorb_field(2,igll,iface) = ty*jacobianw
b_absorb_field(3,igll,iface) = tz*jacobianw
@@ -202,7 +199,100 @@
endif
end subroutine compute_stacey_viscoelastic
+!
+!=====================================================================
+! for elastic solver
+
+! absorbing boundary term for elastic media (Stacey conditions)
+
+ subroutine compute_stacey_viscoelastic_bpwf(NSPEC_AB, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ abs_boundary_ijk,abs_boundary_ispec, &
+ num_abs_boundary_faces, &
+ ispec_is_elastic,SIMULATION_TYPE, &
+ NSTEP,it,NGLOB_ADJOINT,b_accel, &
+ b_num_abs_boundary_faces,b_reclen_field,b_absorb_field)
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+ logical, dimension(NSPEC_AB) :: ispec_is_elastic
+
+! absorbing boundary surface
+ integer :: num_abs_boundary_faces
+ integer :: abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces)
+ integer :: abs_boundary_ispec(num_abs_boundary_faces)
+
+! adjoint simulations
+ integer:: SIMULATION_TYPE
+ integer:: NSTEP,it,NGLOB_ADJOINT
+ integer:: b_num_abs_boundary_faces,b_reclen_field
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLLSQUARE,b_num_abs_boundary_faces):: b_absorb_field
+
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
+
+! local parameters
+ integer :: ispec,iglob,i,j,k,iface,igll
+
+ ! checks if anything to do
+ if( num_abs_boundary_faces == 0 ) return
+
+! adjoint simulations:
+ if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
+ ! reads in absorbing boundary array when first phase is running
+ if( phase_is_inner .eqv. .false. ) then
+ ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newmark scheme
+ ! uses fortran routine
+ !read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
+ !if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
+ ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_field')
+ ! uses c routine for faster reading
+ call read_abs(0,b_absorb_field,b_reclen_field,NSTEP-it+1)
+ endif
+ endif !adjoint
+
+ ! absorbs absorbing-boundary surface using Stacey condition (Clayton & Enquist)
+ do iface=1,num_abs_boundary_faces
+
+ ispec = abs_boundary_ispec(iface)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ if( ispec_is_elastic(ispec) ) then
+
+ ! reference gll points on boundary face
+ do igll = 1,NGLLSQUARE
+
+ ! gets local indices for GLL point
+ i = abs_boundary_ijk(1,igll,iface)
+ j = abs_boundary_ijk(2,igll,iface)
+ k = abs_boundary_ijk(3,igll,iface)
+
+ ! gets velocity
+ iglob=ibool(i,j,k,ispec)
+
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ b_accel(:,iglob) = b_accel(:,iglob) - b_absorb_field(:,igll,iface)
+ endif !adjoint
+
+ enddo
+ endif ! ispec_is_elastic
+ endif ! ispec_is_inner
+ enddo
+
+ end subroutine compute_stacey_viscoelastic_bpwf
+
!---------------------------------------------------------------------------------------
subroutine read_dsm_file(Veloc_dsm_boundary,Tract_dsm_boundary,num_abs_boundary_faces,it_dsm)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90 2013-06-08 12:12:16 UTC (rev 22193)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90 2013-06-09 15:22:30 UTC (rev 22194)
@@ -104,13 +104,30 @@
call it_update_displacement_scheme()
if(.not. GPU_MODE)then
+ if(SIMULATION_TYPE == 3)then
+ if(ELASTIC_SIMULATION .and. ACOUSTIC_SIMULATION)then
+ if( ELASTIC_SIMULATION ) call compute_forces_viscoelastic()
+ if( ACOUSTIC_SIMULATION ) call compute_forces_acoustic()
+ else
+ if( ACOUSTIC_SIMULATION ) call compute_forces_acoustic()
+ if( ELASTIC_SIMULATION ) call compute_forces_viscoelastic()
+ endif
+ else
+ if( ACOUSTIC_SIMULATION ) call compute_forces_acoustic()
+ if( ELASTIC_SIMULATION ) call compute_forces_viscoelastic()
+ endif
+
+ if(SIMULATION_TYPE == 3)then
+ ! acoustic solver
+ ! (needs to be done after elastic one)
+ if( ACOUSTIC_SIMULATION ) call compute_forces_acoustic_bpwf()
+ ! elastic solver
+ ! (needs to be done first, before poroelastic one)
+ if( ELASTIC_SIMULATION ) call compute_forces_viscoelastic_bpwf()
+ endif
+ else
! acoustic solver
- ! (needs to be done first, before elastic one)
- if( ACOUSTIC_SIMULATION ) call compute_forces_acoustic()
- ! elastic solver
- ! (needs to be done first, before poroelastic one)
- if( ELASTIC_SIMULATION ) call compute_forces_viscoelastic()
- else
+ ! (needs to be done after elastic one)
if( ACOUSTIC_SIMULATION ) call compute_forces_acoustic_GPU()
! elastic solver
! (needs to be done first, before poroelastic one)
@@ -917,10 +934,10 @@
if(GPU_MODE) &
call transfer_b_fields_att_to_device(Mesh_pointer, &
b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz,size(b_R_xx), &
-!ZN b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz,size(b_R_xx), & ! please change the above line with this
+!!! b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz,size(b_R_xx), & ! please change the above line with this
b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz, &
-!ZN b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz, &
-!ZN ! please change the above line with this
+!!! b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz, &
+!!! ! please change the above line with this
size(b_epsilondev_xx))
endif
@@ -994,9 +1011,9 @@
! attenuation arrays
call transfer_b_fields_att_to_device(Mesh_pointer, &
b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz,size(b_R_xx), &
-!ZN b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz,size(b_R_xx), &
+!!! b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz,size(b_R_xx), &
b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz, &
-!ZN b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz, &
+!!! b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz, &
size(b_epsilondev_xx))
endif
endif
@@ -1034,9 +1051,9 @@
! attenuation arrays
call transfer_fields_att_from_device(Mesh_pointer, &
R_xx,R_yy,R_xy,R_xz,R_yz,size(R_xx), &
-!ZN R_trace,R_xx,R_yy,R_xy,R_xz,R_yz,size(R_xx), &
+!!! R_trace,R_xx,R_yy,R_xy,R_xz,R_yz,size(R_xx), &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
-!ZN epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
+!!! epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
size(epsilondev_xx))
endif
@@ -1124,9 +1141,9 @@
if (ATTENUATION) &
call transfer_fields_att_from_device(Mesh_pointer, &
R_xx,R_yy,R_xy,R_xz,R_yz,size(R_xx), &
-!ZN R_trace,R_xx,R_yy,R_xy,R_xz,R_yz,size(R_xx), &
+!!! R_trace,R_xx,R_yy,R_xy,R_xz,R_yz,size(R_xx), &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
-!ZN epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
+!!! epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
size(epsilondev_xx))
endif
More information about the CIG-COMMITS
mailing list