[cig-commits] r17141 - seismo/3D/SPECFEM3D/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Sun Aug 29 08:46:16 PDT 2010
Author: danielpeter
Date: 2010-08-29 08:46:15 -0700 (Sun, 29 Aug 2010)
New Revision: 17141
Added:
seismo/3D/SPECFEM3D/trunk/compute_kernels.f90
Modified:
seismo/3D/SPECFEM3D/trunk/Makefile.in
seismo/3D/SPECFEM3D/trunk/combine_surf_data.f90
seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90
seismo/3D/SPECFEM3D/trunk/compute_boundary_kernel.f90
seismo/3D/SPECFEM3D/trunk/compute_forces_elastic.f90
seismo/3D/SPECFEM3D/trunk/compute_forces_elastic_Dev.f90
seismo/3D/SPECFEM3D/trunk/compute_forces_elastic_noDev.f90
seismo/3D/SPECFEM3D/trunk/compute_interpolated_dva.f90
seismo/3D/SPECFEM3D/trunk/constants.h.in
seismo/3D/SPECFEM3D/trunk/finalize_simulation.f90
seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/iterate_time.f90
seismo/3D/SPECFEM3D/trunk/prepare_timerun.f90
seismo/3D/SPECFEM3D/trunk/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/save_adjoint_kernels.f90
seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/write_PNM_GIF_data.f90
Log:
updated elastic kernel calculations; adding new file compute_kernels.f90
Modified: seismo/3D/SPECFEM3D/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/Makefile.in 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/Makefile.in 2010-08-29 15:46:15 UTC (rev 17141)
@@ -138,6 +138,7 @@
$O/specfem3D_par.o \
$O/PML_init.o \
$O/compute_boundary_kernel.o \
+ $O/compute_kernels.o \
$O/compute_forces_acoustic.o \
$O/compute_forces_acoustic_pot.o \
$O/compute_forces_acoustic_PML.o \
@@ -196,6 +197,8 @@
@COND_PYRE_FALSE@ specfem3D \
@COND_PYRE_FALSE@ check_buffers_2D \
@COND_PYRE_FALSE@ combine_AVS_DX \
+ at COND_PYRE_FALSE@ combine_vol_data \
+ at COND_PYRE_FALSE@ combine_surf_data \
@COND_PYRE_FALSE@ convolve_source_timefunction \
@COND_PYRE_FALSE@ $(EMPTY_MACRO)
@@ -573,6 +576,9 @@
$O/compute_boundary_kernel.o: constants.h compute_boundary_kernel.f90
${FCCOMPILE_CHECK} -c -o $O/compute_boundary_kernel.o compute_boundary_kernel.f90
+$O/compute_kernels.o: constants.h compute_kernels.f90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_kernels.o compute_kernels.f90
+
$O/combine_vol_data.o: constants.h combine_vol_data.f90
${FCCOMPILE_CHECK} -c -o $O/combine_vol_data.o combine_vol_data.f90
Modified: seismo/3D/SPECFEM3D/trunk/combine_surf_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/combine_surf_data.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/combine_surf_data.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -39,19 +39,28 @@
integer i,j,k,ispec, ios, it
integer iproc, proc1, proc2, num_node, node_list(300), nspec, nglob
integer np, ne, npp, nee, npoint, nelement, njunk, n1, n2, n3, n4
- integer ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
integer numpoin, iglob1, iglob2, iglob3, iglob4, iglob
- logical mask_ibool(NGLOB_AB)
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: data_3D
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: data_2D
- real(kind=CUSTOM_REAL),dimension(NGLOB_AB) :: xstore, ystore, zstore
+! integer ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
+! logical mask_ibool(NGLOB_AB)
+! real(kind=CUSTOM_REAL),dimension(NGLOB_AB) :: xstore, ystore, zstore
+ ! mesh coordinates
+ real(kind=CUSTOM_REAL),dimension(:),allocatable :: xstore, ystore, zstore
+ integer, dimension(:,:,:,:),allocatable :: ibool
+ logical, dimension(:),allocatable :: mask_ibool
+ integer :: NSPEC_AB,NGLOB_AB
+
real x, y, z
real, dimension(:,:,:,:), allocatable :: dat3D
real, dimension(:,:,:), allocatable :: dat2D
character(len=256) :: sline, arg(8), filename, indir, outdir, prname, surfname
character(len=256) :: mesh_file, local_file, local_data_file, local_ibool_file
character(len=256) :: local_ibool_surf_file
- integer :: num_ibool(NGLOB_AB)
+
+! integer :: num_ibool(NGLOB_AB)
+ integer,dimension(:),allocatable :: num_ibool
+
logical :: HIGH_RESOLUTION_MESH, FILE_ARRAY_IS_3D
integer :: ires, nspec_surf, npoint1, npoint2, ispec_surf, inx, iny, idim
integer,dimension(:), allocatable :: ibelm_surf
@@ -132,8 +141,8 @@
mesh_file = trim(outdir) // '/' // trim(filename)//'.surf'
call open_file(trim(mesh_file)//char(0))
- nspec = NSPEC_AB
- nglob = NGLOB_AB
+! nspec = NSPEC_AB
+! nglob = NGLOB_AB
np = 0
@@ -147,6 +156,21 @@
print *, 'Reading slice ', iproc
write(prname,'(a,i6.6,a)') trim(indir)//'/proc',iproc,'_'
+ ! gets number of elements and global points for this partition
+ open(unit=27,file=prname(1:len_trim(prname))//'external_mesh.bin',&
+ status='old',action='read',form='unformatted',iostat=ios)
+ read(27) NSPEC_AB
+ read(27) NGLOB_AB
+ close(27)
+ nspec = NSPEC_AB
+ nglob = NGLOB_AB
+
+ ! allocates arrays
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(mask_ibool(NGLOB_AB))
+ allocate(num_ibool(NGLOB_AB))
+ allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB))
+
! surface file
local_ibool_surf_file = trim(prname) // 'ibelm_' //trim(surfname)// '.bin'
open(unit = 28,file = trim(local_ibool_surf_file),status='old', iostat = ios, form='unformatted')
@@ -266,6 +290,9 @@
if (numpoin /= npoint) stop 'Error: number of points are not consistent'
np = np + npoint
+ ! frees arrays
+ deallocate(ibool,mask_ibool,num_ibool,xstore,ystore,zstore)
+
enddo ! all slices for points
if (np /= npp) stop 'Error: Number of total points are not consistent'
@@ -282,6 +309,21 @@
print *, 'Reading slice ', iproc
write(prname,'(a,i6.6,a)') trim(indir)//'/proc',iproc,'_'
+ ! gets number of elements and global points for this partition
+ open(unit=27,file=prname(1:len_trim(prname))//'external_mesh.bin',&
+ status='old',action='read',form='unformatted',iostat=ios)
+ read(27) NSPEC_AB
+ read(27) NGLOB_AB
+ close(27)
+ nspec = NSPEC_AB
+ nglob = NGLOB_AB
+
+ ! allocates arrays
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(mask_ibool(NGLOB_AB))
+ allocate(num_ibool(NGLOB_AB))
+
+
np = npoint * (it-1)
! surface file
@@ -351,6 +393,9 @@
enddo
ne = ne + nelement
+ ! frees arrays
+ deallocate(ibool,mask_ibool,num_ibool)
+
enddo ! num_node
if (ne /= nee) stop 'Number of total elements are not consistent'
print *, 'Total number of elements: ', ne
Modified: seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -196,7 +196,8 @@
integer icomp, itime, i, j, k, ios
double precision :: junk
- character(len=3),dimension(NDIM) :: comp = (/ "BHN", "BHE", "BHZ" /)
+ ! note: should have some order as orientation in write_seismograms_to_file()
+ character(len=3),dimension(NDIM) :: comp = (/ "BHE", "BHN", "BHZ" /)
character(len=256) :: filename
!adj_sourcearray(:,:,:,:,:) = 0.
Modified: seismo/3D/SPECFEM3D/trunk/compute_boundary_kernel.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_boundary_kernel.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/compute_boundary_kernel.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -23,9 +23,8 @@
!
!=====================================================================
-subroutine compute_boundary_kernel()
+ subroutine compute_boundary_kernel()
-
! isotropic topography kernel computation
! compare with Tromp et al. (2005), eq. (25), or see Liu & Tromp (2008), eq. (65)
Modified: seismo/3D/SPECFEM3D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_forces_elastic.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/compute_forces_elastic.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -48,32 +48,30 @@
endif
! elastic term
- if(USE_DEVILLE_PRODUCTS) then
+ if(USE_DEVILLE_PRODUCTS) then
call compute_forces_elastic_Dev(iphase, NSPEC_AB,NGLOB_AB,displ,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
ATTENUATION,USE_OLSEN_ATTENUATION, &
- one_minus_sum_beta,factor_common,alphaval,betaval,gammaval, &
- NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ one_minus_sum_beta,factor_common, &
+ alphaval,betaval,gammaval, &
+ NSPEC_ATTENUATION_AB, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy, &
- epsilondev_xz,epsilondev_yz,iflag_attenuation_store, &
+ epsilondev_xz,epsilondev_yz,epsilon_trace_over_3, &
+ iflag_attenuation_store, &
rho_vs,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,NGLOB_ADJOINT,NSPEC_ADJOINT, &
- b_displ,b_accel,kappa_kl,mu_kl,deltat, &
- NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ATT_AND_KERNEL,&
+ SIMULATION_TYPE, COMPUTE_AND_STORE_STRAIN,NSPEC_STRAIN_ONLY, &
+ NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ADJOINT,&
is_moho_top,is_moho_bot, &
- dsdx_top,dsdx_bot,b_dsdx_top,b_dsdx_bot, &
+ dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
- b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
- b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
- b_epsilondev_xz,b_epsilondev_yz, &
- b_alphaval,b_betaval,b_gammaval,&
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic,&
phase_ispec_inner_elastic )
else
@@ -84,29 +82,90 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
ATTENUATION,USE_OLSEN_ATTENUATION,&
- one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
- NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ one_minus_sum_beta,factor_common, &
+ alphaval,betaval,gammaval,&
+ NSPEC_ATTENUATION_AB, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy,&
- epsilondev_xz,epsilondev_yz,iflag_attenuation_store,&
+ epsilondev_xz,epsilondev_yz,epsilon_trace_over_3, &
+ iflag_attenuation_store,&
rho_vs,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,NGLOB_ADJOINT,NSPEC_ADJOINT, &
- b_displ,b_accel,kappa_kl,mu_kl,deltat, &
- NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ATT_AND_KERNEL,&
+ SIMULATION_TYPE,COMPUTE_AND_STORE_STRAIN,NSPEC_STRAIN_ONLY, &
+ NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ADJOINT,&
is_moho_top,is_moho_bot, &
- dsdx_top,dsdx_bot,b_dsdx_top,b_dsdx_bot, &
+ dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
+ num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic,&
+ phase_ispec_inner_elastic )
+ endif
+
+ ! adjoint simulations: backward/reconstructed wavefield
+ if( SIMULATION_TYPE == 3 ) then
+ if(USE_DEVILLE_PRODUCTS) then
+ call compute_forces_elastic_Dev(iphase, NSPEC_AB,NGLOB_AB, &
+ b_displ,b_accel, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ kappastore,mustore,jacobian,ibool, &
+ ATTENUATION,USE_OLSEN_ATTENUATION, &
+ one_minus_sum_beta,factor_common, &
+ b_alphaval,b_betaval,b_gammaval, &
+ NSPEC_ATTENUATION_AB, &
b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
- b_epsilondev_xz,b_epsilondev_yz, &
+ b_epsilondev_xz,b_epsilondev_yz,b_epsilon_trace_over_3, &
+ iflag_attenuation_store, &
+ rho_vs,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 )
+ else
+ call compute_forces_elastic_noDev( iphase, NSPEC_AB,NGLOB_AB,&
+ b_displ,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,USE_OLSEN_ATTENUATION,&
+ one_minus_sum_beta,factor_common, &
b_alphaval,b_betaval,b_gammaval,&
+ NSPEC_ATTENUATION_AB, &
+ b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
+ b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,&
+ b_epsilondev_xz,b_epsilondev_yz,b_epsilon_trace_over_3, &
+ iflag_attenuation_store,&
+ rho_vs,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 )
+
+ endif
endif
+
+
! adds elastic absorbing boundary term to acceleration (Stacey conditions)
if(ABSORBING_CONDITIONS) &
call compute_stacey_elastic(NSPEC_AB,NGLOB_AB,accel, &
Modified: seismo/3D/SPECFEM3D/trunk/compute_forces_elastic_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_forces_elastic_Dev.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/compute_forces_elastic_Dev.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -33,25 +33,21 @@
kappastore,mustore,jacobian,ibool, &
ATTENUATION,USE_OLSEN_ATTENUATION, &
one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
- NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ NSPEC_ATTENUATION_AB, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy, &
- epsilondev_xz,epsilondev_yz,iflag_attenuation_store, &
- rho_vs, &
+ epsilondev_xz,epsilondev_yz,epsilon_trace_over_3, &
+ iflag_attenuation_store,rho_vs, &
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,NGLOB_ADJOINT,NSPEC_ADJOINT, &
- b_displ,b_accel,kappa_kl,mu_kl,deltat, &
- NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ATT_AND_KERNEL, &
+ SIMULATION_TYPE,COMPUTE_AND_STORE_STRAIN,NSPEC_STRAIN_ONLY, &
+ NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ADJOINT, &
is_moho_top,is_moho_bot, &
- dsdx_top,dsdx_bot,b_dsdx_top,b_dsdx_bot, &
+ dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
- b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
- b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
- b_epsilondev_xz,b_epsilondev_yz, &
- b_alphaval,b_betaval,b_gammaval, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic,&
phase_ispec_inner_elastic)
@@ -63,8 +59,6 @@
ONE_THIRD,FOUR_THIRDS,m1,m2
implicit none
- !include "constants.h"
-
integer :: NSPEC_AB,NGLOB_AB
! displacement and acceleration
@@ -83,20 +77,21 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
-! communication overlap
- !logical, dimension(NSPEC_AB) :: ispec_is_inner
- !logical :: phase_is_inner
-
! memory variables and standard linear solids for attenuation
logical :: ATTENUATION,USE_OLSEN_ATTENUATION
+ logical :: COMPUTE_AND_STORE_STRAIN
+ integer :: NSPEC_STRAIN_ONLY, NSPEC_ADJOINT
integer :: NSPEC_ATTENUATION_AB
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: iflag_attenuation_store
real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION) :: one_minus_sum_beta
- real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: factor_common, alphaval,betaval,gammaval
+ real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: factor_common, &
+ alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS) :: &
- R_xx,R_yy,R_xy,R_xz,R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB) :: &
+ R_xx,R_yy,R_xy,R_xz,R_yz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: epsilon_trace_over_3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vs
@@ -109,38 +104,20 @@
c34store,c35store,c36store,c44store,c45store,c46store, &
c55store,c56store,c66store
- !logical,dimension(NSPEC_AB) :: ispec_is_elastic
integer :: iphase
integer :: num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic
integer, dimension(num_phase_ispec_elastic,2) :: phase_ispec_inner_elastic
! adjoint simulations
integer :: SIMULATION_TYPE
- integer :: NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ATT_AND_KERNEL
- integer :: NGLOB_ADJOINT,NSPEC_ADJOINT
+ integer :: NSPEC_BOUN,NSPEC2D_MOHO
! moho kernel
real(kind=CUSTOM_REAL),dimension(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO):: &
- dsdx_top,dsdx_bot,b_dsdx_top,b_dsdx_bot
+ dsdx_top,dsdx_bot
logical,dimension(NSPEC_BOUN) :: is_moho_top,is_moho_bot
integer :: ispec2D_moho_top, ispec2D_moho_bot
-
- ! adjoint memory variables
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS) :: &
- b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL) :: &
- b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz
- real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: b_alphaval,b_betaval,b_gammaval
-
- ! adjoint wavefields
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT):: b_displ,b_accel
- ! adjoint kernels
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: &
- mu_kl, kappa_kl
- real(kind=CUSTOM_REAL) :: deltat
-
-!adjoint
-
+
! local parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc, &
newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3
@@ -184,7 +161,7 @@
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
real(kind=CUSTOM_REAL) factor_loc,alphaval_loc,betaval_loc,gammaval_loc,Sn,Snp1
- real(kind=CUSTOM_REAL) epsilon_trace_over_3
+ real(kind=CUSTOM_REAL) templ
real(kind=CUSTOM_REAL) vs_val
real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
@@ -205,71 +182,21 @@
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
integer i_SLS,iselected
-
integer ispec,iglob,ispec_p,num_elements
integer i,j,k
- ! adjoint backward arrays
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: b_dummyx_loc,b_dummyy_loc,b_dummyz_loc, &
- b_newtempx1,b_newtempx2,b_newtempx3,b_newtempy1,b_newtempy2,b_newtempy3,b_newtempz1,b_newtempz2,b_newtempz3
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
- b_tempx1,b_tempx2,b_tempx3,b_tempy1,b_tempy2,b_tempy3,b_tempz1,b_tempz2,b_tempz3
- ! backward arrays: manually inline the calls to the Deville et al. (2002) routines
- real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: b_B1_m1_m2_5points,b_B2_m1_m2_5points,b_B3_m1_m2_5points
- real(kind=CUSTOM_REAL), dimension(m1,m2) :: b_C1_m1_m2_5points,b_C2_m1_m2_5points,b_C3_m1_m2_5points
- real(kind=CUSTOM_REAL), dimension(m1,m2) :: b_E1_m1_m2_5points,b_E2_m1_m2_5points,b_E3_m1_m2_5points
- equivalence(b_dummyx_loc,b_B1_m1_m2_5points)
- equivalence(b_dummyy_loc,b_B2_m1_m2_5points)
- equivalence(b_dummyz_loc,b_B3_m1_m2_5points)
- equivalence(b_tempx1,b_C1_m1_m2_5points)
- equivalence(b_tempy1,b_C2_m1_m2_5points)
- equivalence(b_tempz1,b_C3_m1_m2_5points)
- equivalence(b_newtempx1,b_E1_m1_m2_5points)
- equivalence(b_newtempy1,b_E2_m1_m2_5points)
- equivalence(b_newtempz1,b_E3_m1_m2_5points)
- real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: &
- b_A1_mxm_m2_m1_5points,b_A2_mxm_m2_m1_5points,b_A3_mxm_m2_m1_5points
- real(kind=CUSTOM_REAL), dimension(m2,m1) :: &
- b_C1_mxm_m2_m1_5points,b_C2_mxm_m2_m1_5points,b_C3_mxm_m2_m1_5points
- real(kind=CUSTOM_REAL), dimension(m2,m1) :: &
- b_E1_mxm_m2_m1_5points,b_E2_mxm_m2_m1_5points,b_E3_mxm_m2_m1_5points
- equivalence(b_dummyx_loc,b_A1_mxm_m2_m1_5points)
- equivalence(b_dummyy_loc,b_A2_mxm_m2_m1_5points)
- equivalence(b_dummyz_loc,b_A3_mxm_m2_m1_5points)
- equivalence(b_tempx3,b_C1_mxm_m2_m1_5points)
- equivalence(b_tempy3,b_C2_mxm_m2_m1_5points)
- equivalence(b_tempz3,b_C3_mxm_m2_m1_5points)
- equivalence(b_newtempx3,b_E1_mxm_m2_m1_5points)
- equivalence(b_newtempy3,b_E2_mxm_m2_m1_5points)
- equivalence(b_newtempz3,b_E3_mxm_m2_m1_5points)
- real(kind=CUSTOM_REAL):: dsxx,dsxy,dsxz,dsyy,dsyz,dszz
- real(kind=CUSTOM_REAL):: b_duxdxl,b_duxdyl,b_duxdzl,b_duydxl,b_duydyl,b_duydzl,b_duzdxl,b_duzdyl,b_duzdzl
- real(kind=CUSTOM_REAL):: b_duxdxl_plus_duydyl,b_duxdxl_plus_duzdzl,b_duydyl_plus_duzdzl
- real(kind=CUSTOM_REAL):: b_duxdyl_plus_duydxl,b_duzdxl_plus_duxdzl,b_duzdyl_plus_duydzl
- real(kind=CUSTOM_REAL):: b_dsxx,b_dsxy,b_dsxz,b_dsyy,b_dsyz,b_dszz
- real(kind=CUSTOM_REAL):: b_sigma_xx,b_sigma_yy,b_sigma_zz,b_sigma_xy,b_sigma_xz,b_sigma_yz
- real(kind=CUSTOM_REAL):: kappa_k, mu_k
- ! local adjoint attenuation
- real(kind=CUSTOM_REAL) b_alphaval_loc,b_betaval_loc,b_gammaval_loc,b_Sn,b_Snp1
- real(kind=CUSTOM_REAL) b_epsilon_trace_over_3
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: b_epsilondev_xx_loc, &
- b_epsilondev_yy_loc, b_epsilondev_xy_loc, b_epsilondev_xz_loc, b_epsilondev_yz_loc
- real(kind=CUSTOM_REAL) b_R_xx_val,b_R_yy_val
- ! adjoint
+! real(kind=CUSTOM_REAL):: dsxx,dsxy,dsxz,dsyy,dsyz,dszz
+ ! choses inner/outer elements
if( iphase == 1 ) then
num_elements = nspec_outer_elastic
else
num_elements = nspec_inner_elastic
endif
-! loops over all elements
-! do ispec = 1,NSPEC_AB
-! if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
-! if( ispec_is_elastic(ispec) ) then
-
do ispec_p = 1,num_elements
+ ! returns element id from stored element list
ispec = phase_ispec_inner_elastic(ispec_p,iphase)
! adjoint simulations: moho kernel
@@ -289,13 +216,6 @@
dummyx_loc(i,j,k) = displ(1,iglob)
dummyy_loc(i,j,k) = displ(2,iglob)
dummyz_loc(i,j,k) = displ(3,iglob)
-
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- b_dummyx_loc(i,j,k) = b_displ(1,iglob)
- b_dummyy_loc(i,j,k) = b_displ(2,iglob)
- b_dummyz_loc(i,j,k) = b_displ(3,iglob)
- endif
enddo
enddo
enddo
@@ -321,26 +241,6 @@
hprime_xx(i,3)*B3_m1_m2_5points(3,j) + &
hprime_xx(i,4)*B3_m1_m2_5points(4,j) + &
hprime_xx(i,5)*B3_m1_m2_5points(5,j)
-
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- b_C1_m1_m2_5points(i,j) = hprime_xx(i,1)*b_B1_m1_m2_5points(1,j) + &
- hprime_xx(i,2)*b_B1_m1_m2_5points(2,j) + &
- hprime_xx(i,3)*b_B1_m1_m2_5points(3,j) + &
- hprime_xx(i,4)*b_B1_m1_m2_5points(4,j) + &
- hprime_xx(i,5)*b_B1_m1_m2_5points(5,j)
- b_C2_m1_m2_5points(i,j) = hprime_xx(i,1)*b_B2_m1_m2_5points(1,j) + &
- hprime_xx(i,2)*b_B2_m1_m2_5points(2,j) + &
- hprime_xx(i,3)*b_B2_m1_m2_5points(3,j) + &
- hprime_xx(i,4)*b_B2_m1_m2_5points(4,j) + &
- hprime_xx(i,5)*b_B2_m1_m2_5points(5,j)
- b_C3_m1_m2_5points(i,j) = hprime_xx(i,1)*b_B3_m1_m2_5points(1,j) + &
- hprime_xx(i,2)*b_B3_m1_m2_5points(2,j) + &
- hprime_xx(i,3)*b_B3_m1_m2_5points(3,j) + &
- hprime_xx(i,4)*b_B3_m1_m2_5points(4,j) + &
- hprime_xx(i,5)*b_B3_m1_m2_5points(5,j)
- endif ! adjoint
-
enddo
enddo
@@ -365,25 +265,6 @@
dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
dummyz_loc(i,5,k)*hprime_xxT(5,j)
-
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- b_tempx2(i,j,k) = b_dummyx_loc(i,1,k)*hprime_xxT(1,j) + &
- b_dummyx_loc(i,2,k)*hprime_xxT(2,j) + &
- b_dummyx_loc(i,3,k)*hprime_xxT(3,j) + &
- b_dummyx_loc(i,4,k)*hprime_xxT(4,j) + &
- b_dummyx_loc(i,5,k)*hprime_xxT(5,j)
- b_tempy2(i,j,k) = b_dummyy_loc(i,1,k)*hprime_xxT(1,j) + &
- b_dummyy_loc(i,2,k)*hprime_xxT(2,j) + &
- b_dummyy_loc(i,3,k)*hprime_xxT(3,j) + &
- b_dummyy_loc(i,4,k)*hprime_xxT(4,j) + &
- b_dummyy_loc(i,5,k)*hprime_xxT(5,j)
- b_tempz2(i,j,k) = b_dummyz_loc(i,1,k)*hprime_xxT(1,j) + &
- b_dummyz_loc(i,2,k)*hprime_xxT(2,j) + &
- b_dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
- b_dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
- b_dummyz_loc(i,5,k)*hprime_xxT(5,j)
- endif ! adjoint
enddo
enddo
enddo
@@ -406,25 +287,6 @@
A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
-
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- b_C1_mxm_m2_m1_5points(i,j) = b_A1_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
- b_A1_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
- b_A1_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
- b_A1_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
- b_A1_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
- b_C2_mxm_m2_m1_5points(i,j) = b_A2_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
- b_A2_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
- b_A2_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
- b_A2_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
- b_A2_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
- b_C3_mxm_m2_m1_5points(i,j) = b_A3_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
- b_A3_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
- b_A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
- b_A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
- b_A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
- endif ! adjoint
enddo
enddo
@@ -455,6 +317,31 @@
duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+ ! save strain on the Moho boundary
+ if (SAVE_MOHO_MESH ) then
+ if (is_moho_top(ispec)) then
+ dsdx_top(1,1,i,j,k,ispec2D_moho_top) = duxdxl
+ dsdx_top(1,2,i,j,k,ispec2D_moho_top) = duxdyl
+ dsdx_top(1,3,i,j,k,ispec2D_moho_top) = duxdzl
+ dsdx_top(2,1,i,j,k,ispec2D_moho_top) = duydxl
+ dsdx_top(2,2,i,j,k,ispec2D_moho_top) = duydyl
+ dsdx_top(2,3,i,j,k,ispec2D_moho_top) = duydzl
+ dsdx_top(3,1,i,j,k,ispec2D_moho_top) = duzdxl
+ dsdx_top(3,2,i,j,k,ispec2D_moho_top) = duzdyl
+ dsdx_top(3,3,i,j,k,ispec2D_moho_top) = duzdzl
+ else if (is_moho_bot(ispec)) then
+ dsdx_bot(1,1,i,j,k,ispec2D_moho_bot) = duxdxl
+ dsdx_bot(1,2,i,j,k,ispec2D_moho_bot) = duxdyl
+ dsdx_bot(1,3,i,j,k,ispec2D_moho_bot) = duxdzl
+ dsdx_bot(2,1,i,j,k,ispec2D_moho_bot) = duydxl
+ dsdx_bot(2,2,i,j,k,ispec2D_moho_bot) = duydyl
+ dsdx_bot(2,3,i,j,k,ispec2D_moho_bot) = duydzl
+ dsdx_bot(3,1,i,j,k,ispec2D_moho_bot) = duzdxl
+ dsdx_bot(3,2,i,j,k,ispec2D_moho_bot) = duzdyl
+ dsdx_bot(3,3,i,j,k,ispec2D_moho_bot) = duzdzl
+ endif
+ endif
+
! precompute some sums to save CPU time
duxdxl_plus_duydyl = duxdxl + duydyl
duxdxl_plus_duzdzl = duxdxl + duzdzl
@@ -463,122 +350,22 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
+
kappal = kappastore(i,j,k,ispec)
mul = mustore(i,j,k,ispec)
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- ! save strain on the Moho boundary
- if (SAVE_MOHO_MESH ) then
- if (is_moho_top(ispec)) then
- dsdx_top(1,1,i,j,k,ispec2D_moho_top) = duxdxl
- dsdx_top(1,2,i,j,k,ispec2D_moho_top) = duxdyl
- dsdx_top(1,3,i,j,k,ispec2D_moho_top) = duxdzl
- dsdx_top(2,1,i,j,k,ispec2D_moho_top) = duydxl
- dsdx_top(2,2,i,j,k,ispec2D_moho_top) = duydyl
- dsdx_top(2,3,i,j,k,ispec2D_moho_top) = duydzl
- dsdx_top(3,1,i,j,k,ispec2D_moho_top) = duzdxl
- dsdx_top(3,2,i,j,k,ispec2D_moho_top) = duzdyl
- dsdx_top(3,3,i,j,k,ispec2D_moho_top) = duzdzl
- else if (is_moho_bot(ispec)) then
- dsdx_bot(1,1,i,j,k,ispec2D_moho_bot) = duxdxl
- dsdx_bot(1,2,i,j,k,ispec2D_moho_bot) = duxdyl
- dsdx_bot(1,3,i,j,k,ispec2D_moho_bot) = duxdzl
- dsdx_bot(2,1,i,j,k,ispec2D_moho_bot) = duydxl
- dsdx_bot(2,2,i,j,k,ispec2D_moho_bot) = duydyl
- dsdx_bot(2,3,i,j,k,ispec2D_moho_bot) = duydzl
- dsdx_bot(3,1,i,j,k,ispec2D_moho_bot) = duzdxl
- dsdx_bot(3,2,i,j,k,ispec2D_moho_bot) = duzdyl
- dsdx_bot(3,3,i,j,k,ispec2D_moho_bot) = duzdzl
- endif
- endif
-
- dsxx = duxdxl
- dsxy = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
- dsxz = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
- dsyy = duydyl
- dsyz = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
- dszz = duzdzl
-
- b_duxdxl = xixl*b_tempx1(i,j,k) + etaxl*b_tempx2(i,j,k) + gammaxl*b_tempx3(i,j,k)
- b_duxdyl = xiyl*b_tempx1(i,j,k) + etayl*b_tempx2(i,j,k) + gammayl*b_tempx3(i,j,k)
- b_duxdzl = xizl*b_tempx1(i,j,k) + etazl*b_tempx2(i,j,k) + gammazl*b_tempx3(i,j,k)
- b_duydxl = xixl*b_tempy1(i,j,k) + etaxl*b_tempy2(i,j,k) + gammaxl*b_tempy3(i,j,k)
- b_duydyl = xiyl*b_tempy1(i,j,k) + etayl*b_tempy2(i,j,k) + gammayl*b_tempy3(i,j,k)
- b_duydzl = xizl*b_tempy1(i,j,k) + etazl*b_tempy2(i,j,k) + gammazl*b_tempy3(i,j,k)
- b_duzdxl = xixl*b_tempz1(i,j,k) + etaxl*b_tempz2(i,j,k) + gammaxl*b_tempz3(i,j,k)
- b_duzdyl = xiyl*b_tempz1(i,j,k) + etayl*b_tempz2(i,j,k) + gammayl*b_tempz3(i,j,k)
- b_duzdzl = xizl*b_tempz1(i,j,k) + etazl*b_tempz2(i,j,k) + gammazl*b_tempz3(i,j,k)
-
- b_duxdxl_plus_duydyl = b_duxdxl + b_duydyl
- b_duxdxl_plus_duzdzl = b_duxdxl + b_duzdzl
- b_duydyl_plus_duzdzl = b_duydyl + b_duzdzl
- b_duxdyl_plus_duydxl = b_duxdyl + b_duydxl
- b_duzdxl_plus_duxdzl = b_duzdxl + b_duxdzl
- b_duzdyl_plus_duydzl = b_duzdyl + b_duydzl
-
- b_dsxx = b_duxdxl
- b_dsxy = 0.5_CUSTOM_REAL * b_duxdyl_plus_duydxl
- b_dsxz = 0.5_CUSTOM_REAL * b_duzdxl_plus_duxdzl
- b_dsyy = b_duydyl
- b_dsyz = 0.5_CUSTOM_REAL * b_duzdyl_plus_duydzl
- b_dszz = b_duzdzl
-
- ! isotropic adjoint kernels: bulk (kappa) and shear (mu) kernels
- kappa_k = (duxdxl + duydyl + duzdzl) * (b_duxdxl + b_duydyl + b_duzdzl)
- mu_k = dsxx * b_dsxx + dsyy * b_dsyy + dszz * b_dszz + &
- 2._CUSTOM_REAL * (dsxy * b_dsxy + dsxz * b_dsxz + dsyz * b_dsyz) &
- - ONE_THIRD * kappa_k
-
- kappa_kl(i,j,k,ispec) = kappa_kl(i,j,k,ispec) + deltat * kappa_k
- mu_kl(i,j,k,ispec) = mu_kl(i,j,k,ispec) + 2._CUSTOM_REAL * deltat * mu_k
-
- if (SAVE_MOHO_MESH) then
- if (is_moho_top(ispec)) then
- b_dsdx_top(1,1,i,j,k,ispec2D_moho_top) = b_duxdxl
- b_dsdx_top(1,2,i,j,k,ispec2D_moho_top) = b_duxdyl
- b_dsdx_top(1,3,i,j,k,ispec2D_moho_top) = b_duxdzl
- b_dsdx_top(2,1,i,j,k,ispec2D_moho_top) = b_duydxl
- b_dsdx_top(2,2,i,j,k,ispec2D_moho_top) = b_duydyl
- b_dsdx_top(2,3,i,j,k,ispec2D_moho_top) = b_duydzl
- b_dsdx_top(3,1,i,j,k,ispec2D_moho_top) = b_duzdxl
- b_dsdx_top(3,2,i,j,k,ispec2D_moho_top) = b_duzdyl
- b_dsdx_top(3,3,i,j,k,ispec2D_moho_top) = b_duzdzl
- else if (is_moho_bot(ispec)) then
- b_dsdx_bot(1,1,i,j,k,ispec2D_moho_bot) = b_duxdxl
- b_dsdx_bot(1,2,i,j,k,ispec2D_moho_bot) = b_duxdyl
- b_dsdx_bot(1,3,i,j,k,ispec2D_moho_bot) = b_duxdzl
- b_dsdx_bot(2,1,i,j,k,ispec2D_moho_bot) = b_duydxl
- b_dsdx_bot(2,2,i,j,k,ispec2D_moho_bot) = b_duydyl
- b_dsdx_bot(2,3,i,j,k,ispec2D_moho_bot) = b_duydzl
- b_dsdx_bot(3,1,i,j,k,ispec2D_moho_bot) = b_duzdxl
- b_dsdx_bot(3,2,i,j,k,ispec2D_moho_bot) = b_duzdyl
- b_dsdx_bot(3,3,i,j,k,ispec2D_moho_bot) = b_duzdzl
- endif
- endif
- endif ! adjoint
-
-
! attenuation
if(ATTENUATION) then
- ! compute deviatoric strain
- epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_xx_loc(i,j,k) = duxdxl - epsilon_trace_over_3
- epsilondev_yy_loc(i,j,k) = duydyl - epsilon_trace_over_3
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
-
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_epsilon_trace_over_3 = ONE_THIRD * (b_duxdxl + b_duydyl + b_duzdzl)
- b_epsilondev_xx_loc(i,j,k) = b_duxdxl - b_epsilon_trace_over_3
- b_epsilondev_yy_loc(i,j,k) = b_duydyl - b_epsilon_trace_over_3
- b_epsilondev_xy_loc(i,j,k) = 0.5 * b_duxdyl_plus_duydxl
- b_epsilondev_xz_loc(i,j,k) = 0.5 * b_duzdxl_plus_duxdzl
- b_epsilondev_yz_loc(i,j,k) = 0.5 * b_duzdyl_plus_duydzl
- endif ! adjoint
-
! uses scaling rule similar to Olsen et al. (2003) or mesh flag
if(USE_OLSEN_ATTENUATION) then
vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
@@ -616,18 +403,6 @@
c55 = c55store(i,j,k,ispec)
c56 = c56store(i,j,k,ispec)
c66 = c66store(i,j,k,ispec)
- !if(ATTENUATION .and. not_fully_in_bedrock(ispec)) then
- ! mul = c44
- ! c11 = c11 + FOUR_THIRDS * minus_sum_beta * mul
- ! c12 = c12 - TWO_THIRDS * minus_sum_beta * mul
- ! c13 = c13 - TWO_THIRDS * minus_sum_beta * mul
- ! c22 = c22 + FOUR_THIRDS * minus_sum_beta * mul
- ! c23 = c23 - TWO_THIRDS * minus_sum_beta * mul
- ! c33 = c33 + FOUR_THIRDS * minus_sum_beta * mul
- ! c44 = c44 + minus_sum_beta * mul
- ! c55 = c55 + minus_sum_beta * mul
- ! c66 = c66 + minus_sum_beta * mul
- !endif
sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
@@ -642,21 +417,6 @@
sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_sigma_xx = c11*b_duxdxl + c16*b_duxdyl_plus_duydxl + c12*b_duydyl + &
- c15*b_duzdxl_plus_duxdzl + c14*b_duzdyl_plus_duydzl + c13*b_duzdzl
- b_sigma_yy = c12*b_duxdxl + c26*b_duxdyl_plus_duydxl + c22*b_duydyl + &
- c25*b_duzdxl_plus_duxdzl + c24*b_duzdyl_plus_duydzl + c23*b_duzdzl
- b_sigma_zz = c13*b_duxdxl + c36*b_duxdyl_plus_duydxl + c23*b_duydyl + &
- c35*b_duzdxl_plus_duxdzl + c34*b_duzdyl_plus_duydzl + c33*b_duzdzl
- b_sigma_xy = c16*b_duxdxl + c66*b_duxdyl_plus_duydxl + c26*b_duydyl + &
- c56*b_duzdxl_plus_duxdzl + c46*b_duzdyl_plus_duydzl + c36*b_duzdzl
- b_sigma_xz = c15*b_duxdxl + c56*b_duxdyl_plus_duydxl + c25*b_duydyl + &
- c55*b_duzdxl_plus_duxdzl + c45*b_duzdyl_plus_duydzl + c35*b_duzdzl
- b_sigma_yz = c14*b_duxdxl + c46*b_duxdyl_plus_duydxl + c24*b_duydyl + &
- c45*b_duzdxl_plus_duxdzl + c44*b_duzdyl_plus_duydzl + c34*b_duzdzl
- endif ! adjoint
else
! isotropic case
@@ -672,16 +432,6 @@
sigma_xz = mul*duzdxl_plus_duxdzl
sigma_yz = mul*duzdyl_plus_duydzl
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_sigma_xx = lambdalplus2mul*b_duxdxl + lambdal*b_duydyl_plus_duzdzl
- b_sigma_yy = lambdalplus2mul*b_duydyl + lambdal*b_duxdxl_plus_duzdzl
- b_sigma_zz = lambdalplus2mul*b_duzdzl + lambdal*b_duxdxl_plus_duydyl
- b_sigma_xy = mul*b_duxdyl_plus_duydxl
- b_sigma_xz = mul*b_duzdxl_plus_duxdzl
- b_sigma_yz = mul*b_duzdyl_plus_duydzl
- endif !adjoint
-
endif ! ANISOTROPY
! subtract memory variables if attenuation
@@ -695,18 +445,6 @@
sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
-
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_R_xx_val = b_R_xx(i,j,k,ispec,i_sls)
- b_R_yy_val = b_R_yy(i,j,k,ispec,i_sls)
- b_sigma_xx = b_sigma_xx - b_R_xx_val
- b_sigma_yy = b_sigma_yy - b_R_yy_val
- b_sigma_zz = b_sigma_zz + b_R_xx_val + b_R_yy_val
- b_sigma_xy = b_sigma_xy - b_R_xy(i,j,k,ispec,i_sls)
- b_sigma_xz = b_sigma_xz - b_R_xz(i,j,k,ispec,i_sls)
- b_sigma_yz = b_sigma_yz - b_R_yz(i,j,k,ispec,i_sls)
- endif !adjoint
enddo
endif
@@ -723,19 +461,6 @@
tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx1(i,j,k) = jacobianl * (b_sigma_xx*xixl + b_sigma_xy*xiyl + b_sigma_xz*xizl)
- b_tempy1(i,j,k) = jacobianl * (b_sigma_xy*xixl + b_sigma_yy*xiyl + b_sigma_yz*xizl)
- b_tempz1(i,j,k) = jacobianl * (b_sigma_xz*xixl + b_sigma_yz*xiyl + b_sigma_zz*xizl)
- b_tempx2(i,j,k) = jacobianl * (b_sigma_xx*etaxl + b_sigma_xy*etayl + b_sigma_xz*etazl)
- b_tempy2(i,j,k) = jacobianl * (b_sigma_xy*etaxl + b_sigma_yy*etayl + b_sigma_yz*etazl)
- b_tempz2(i,j,k) = jacobianl * (b_sigma_xz*etaxl + b_sigma_yz*etayl + b_sigma_zz*etazl)
- b_tempx3(i,j,k) = jacobianl * (b_sigma_xx*gammaxl + b_sigma_xy*gammayl + b_sigma_xz*gammazl)
- b_tempy3(i,j,k) = jacobianl * (b_sigma_xy*gammaxl + b_sigma_yy*gammayl + b_sigma_yz*gammazl)
- b_tempz3(i,j,k) = jacobianl * (b_sigma_xz*gammaxl + b_sigma_yz*gammayl + b_sigma_zz*gammazl)
- endif !adjoint
-
enddo
enddo
enddo
@@ -761,25 +486,6 @@
hprimewgll_xxT(i,3)*C3_m1_m2_5points(3,j) + &
hprimewgll_xxT(i,4)*C3_m1_m2_5points(4,j) + &
hprimewgll_xxT(i,5)*C3_m1_m2_5points(5,j)
-
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- b_E1_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*b_C1_m1_m2_5points(1,j) + &
- hprimewgll_xxT(i,2)*b_C1_m1_m2_5points(2,j) + &
- hprimewgll_xxT(i,3)*b_C1_m1_m2_5points(3,j) + &
- hprimewgll_xxT(i,4)*b_C1_m1_m2_5points(4,j) + &
- hprimewgll_xxT(i,5)*b_C1_m1_m2_5points(5,j)
- b_E2_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*b_C2_m1_m2_5points(1,j) + &
- hprimewgll_xxT(i,2)*b_C2_m1_m2_5points(2,j) + &
- hprimewgll_xxT(i,3)*b_C2_m1_m2_5points(3,j) + &
- hprimewgll_xxT(i,4)*b_C2_m1_m2_5points(4,j) + &
- hprimewgll_xxT(i,5)*b_C2_m1_m2_5points(5,j)
- b_E3_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*b_C3_m1_m2_5points(1,j) + &
- hprimewgll_xxT(i,2)*b_C3_m1_m2_5points(2,j) + &
- hprimewgll_xxT(i,3)*b_C3_m1_m2_5points(3,j) + &
- hprimewgll_xxT(i,4)*b_C3_m1_m2_5points(4,j) + &
- hprimewgll_xxT(i,5)*b_C3_m1_m2_5points(5,j)
- endif !adjoint
enddo
enddo
@@ -804,25 +510,6 @@
tempz2(i,3,k)*hprimewgll_xx(3,j) + &
tempz2(i,4,k)*hprimewgll_xx(4,j) + &
tempz2(i,5,k)*hprimewgll_xx(5,j)
-
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- b_newtempx2(i,j,k) = b_tempx2(i,1,k)*hprimewgll_xx(1,j) + &
- b_tempx2(i,2,k)*hprimewgll_xx(2,j) + &
- b_tempx2(i,3,k)*hprimewgll_xx(3,j) + &
- b_tempx2(i,4,k)*hprimewgll_xx(4,j) + &
- b_tempx2(i,5,k)*hprimewgll_xx(5,j)
- b_newtempy2(i,j,k) = b_tempy2(i,1,k)*hprimewgll_xx(1,j) + &
- b_tempy2(i,2,k)*hprimewgll_xx(2,j) + &
- b_tempy2(i,3,k)*hprimewgll_xx(3,j) + &
- b_tempy2(i,4,k)*hprimewgll_xx(4,j) + &
- b_tempy2(i,5,k)*hprimewgll_xx(5,j)
- b_newtempz2(i,j,k) = b_tempz2(i,1,k)*hprimewgll_xx(1,j) + &
- b_tempz2(i,2,k)*hprimewgll_xx(2,j) + &
- b_tempz2(i,3,k)*hprimewgll_xx(3,j) + &
- b_tempz2(i,4,k)*hprimewgll_xx(4,j) + &
- b_tempz2(i,5,k)*hprimewgll_xx(5,j)
- endif !adjoint
enddo
enddo
enddo
@@ -845,25 +532,6 @@
C3_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
C3_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
C3_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
-
- ! adjoint simulations
- if( SIMULATION_TYPE == 3 ) then
- b_E1_mxm_m2_m1_5points(i,j) = b_C1_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
- b_C1_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
- b_C1_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
- b_C1_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
- b_C1_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
- b_E2_mxm_m2_m1_5points(i,j) = b_C2_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
- b_C2_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
- b_C2_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
- b_C2_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
- b_C2_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
- b_E3_mxm_m2_m1_5points(i,j) = b_C3_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
- b_C3_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
- b_C3_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
- b_C3_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
- b_C3_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
- endif !adjoint
enddo
enddo
@@ -884,16 +552,6 @@
accel(3,iglob) = accel(3,iglob) - fac1*newtempz1(i,j,k) - &
fac2*newtempz2(i,j,k) - fac3*newtempz3(i,j,k)
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_accel(1,iglob) = b_accel(1,iglob) - fac1*b_newtempx1(i,j,k) - &
- fac2*b_newtempx2(i,j,k) - fac3*b_newtempx3(i,j,k)
- b_accel(2,iglob) = b_accel(2,iglob) - fac1*b_newtempy1(i,j,k) - &
- fac2*b_newtempy2(i,j,k) - fac3*b_newtempy3(i,j,k)
- b_accel(3,iglob) = b_accel(3,iglob) - fac1*b_newtempz1(i,j,k) - &
- fac2*b_newtempz2(i,j,k) - fac3*b_newtempz3(i,j,k)
- endif !adjoint
-
! update memory variables based upon the Runge-Kutta scheme
if(ATTENUATION) then
@@ -941,39 +599,6 @@
R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + &
betaval_loc * Sn + gammaval_loc * Snp1
- !adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_alphaval_loc = b_alphaval(iselected,i_sls)
- b_betaval_loc = b_betaval(iselected,i_sls)
- b_gammaval_loc = b_gammaval(iselected,i_sls)
- ! term in xx
- b_Sn = factor_loc * b_epsilondev_xx(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_xx_loc(i,j,k)
- b_R_xx(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xx(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- ! term in yy
- b_Sn = factor_loc * b_epsilondev_yy(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_yy_loc(i,j,k)
- b_R_yy(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_yy(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- ! term in zz not computed since zero trace
- ! term in xy
- b_Sn = factor_loc * b_epsilondev_xy(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_xy_loc(i,j,k)
- b_R_xy(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xy(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- ! term in xz
- b_Sn = factor_loc * b_epsilondev_xz(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_xz_loc(i,j,k)
- b_R_xz(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xz(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- ! term in yz
- b_Sn = factor_loc * b_epsilondev_yz(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_yz_loc(i,j,k)
- b_R_yz(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_yz(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- endif !adjoint
-
enddo ! end of loop on memory variables
endif ! end attenuation
@@ -983,152 +608,15 @@
enddo
! save deviatoric strain for Runge-Kutta scheme
- if(ATTENUATION) then
+ if ( COMPUTE_AND_STORE_STRAIN ) then
epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_epsilondev_xx(:,:,:,ispec) = b_epsilondev_xx_loc(:,:,:)
- b_epsilondev_yy(:,:,:,ispec) = b_epsilondev_yy_loc(:,:,:)
- b_epsilondev_xy(:,:,:,ispec) = b_epsilondev_xy_loc(:,:,:)
- b_epsilondev_xz(:,:,:,ispec) = b_epsilondev_xz_loc(:,:,:)
- b_epsilondev_yz(:,:,:,ispec) = b_epsilondev_yz_loc(:,:,:)
- endif !adjoint
endif
-! endif ! ispec_is_elastic
-
-! endif ! if (ispec_is_inner(ispec) .eqv. phase_is_inner)
-
enddo ! spectral element loop
end subroutine compute_forces_elastic_Dev
-!
-!-------------------------------------------------------------------------------------------------
-!
-!
-!! subroutines adapted from Deville, Fischer and Mund, High-order methods
-!! for incompressible fluid flow, Cambridge University Press (2002),
-!! pages 386 and 389 and Figure 8.3.1
-!
-! subroutine old_mxm_m1_m2_5points(A,B1,B2,B3,C1,C2,C3)
-!
-! implicit none
-!
-! include "constants.h"
-!
-! real(kind=4), dimension(m1,NGLLX) :: A
-! real(kind=4), dimension(NGLLX,m2) :: B1,B2,B3
-! real(kind=4), dimension(m1,m2) :: C1,C2,C3
-!
-! integer :: i,j
-!
-! do j=1,m2
-! do i=1,m1
-!
-! C1(i,j) = A(i,1)*B1(1,j) + &
-! A(i,2)*B1(2,j) + &
-! A(i,3)*B1(3,j) + &
-! A(i,4)*B1(4,j) + &
-! A(i,5)*B1(5,j)
-!
-! C2(i,j) = A(i,1)*B2(1,j) + &
-! A(i,2)*B2(2,j) + &
-! A(i,3)*B2(3,j) + &
-! A(i,4)*B2(4,j) + &
-! A(i,5)*B2(5,j)
-!
-! C3(i,j) = A(i,1)*B3(1,j) + &
-! A(i,2)*B3(2,j) + &
-! A(i,3)*B3(3,j) + &
-! A(i,4)*B3(4,j) + &
-! A(i,5)*B3(5,j)
-!
-! enddo
-! enddo
-!
-! end subroutine old_mxm_m1_m2_5points
-!
-!!---------
-!
-! subroutine old_mxm_m1_m1_5points(A1,A2,A3,B,C1,C2,C3)
-!
-! implicit none
-!
-! include "constants.h"
-!
-! real(kind=4), dimension(m1,NGLLX) :: A1,A2,A3
-! real(kind=4), dimension(NGLLX,m1) :: B
-! real(kind=4), dimension(m1,m1) :: C1,C2,C3
-!
-! integer :: i,j
-!
-! do j=1,m1
-! do i=1,m1
-!
-! C1(i,j) = A1(i,1)*B(1,j) + &
-! A1(i,2)*B(2,j) + &
-! A1(i,3)*B(3,j) + &
-! A1(i,4)*B(4,j) + &
-! A1(i,5)*B(5,j)
-!
-! C2(i,j) = A2(i,1)*B(1,j) + &
-! A2(i,2)*B(2,j) + &
-! A2(i,3)*B(3,j) + &
-! A2(i,4)*B(4,j) + &
-! A2(i,5)*B(5,j)
-!
-! C3(i,j) = A3(i,1)*B(1,j) + &
-! A3(i,2)*B(2,j) + &
-! A3(i,3)*B(3,j) + &
-! A3(i,4)*B(4,j) + &
-! A3(i,5)*B(5,j)
-!
-! enddo
-! enddo
-!
-! end subroutine old_mxm_m1_m1_5points
-!
-!!---------
-!
-! subroutine old_mxm_m2_m1_5points(A1,A2,A3,B,C1,C2,C3)
-!
-! implicit none
-!
-! include "constants.h"
-!
-! real(kind=4), dimension(m2,NGLLX) :: A1,A2,A3
-! real(kind=4), dimension(NGLLX,m1) :: B
-! real(kind=4), dimension(m2,m1) :: C1,C2,C3
-!
-! integer :: i,j
-!
-! do j=1,m1
-! do i=1,m2
-!
-! C1(i,j) = A1(i,1)*B(1,j) + &
-! A1(i,2)*B(2,j) + &
-! A1(i,3)*B(3,j) + &
-! A1(i,4)*B(4,j) + &
-! A1(i,5)*B(5,j)
-!
-! C2(i,j) = A2(i,1)*B(1,j) + &
-! A2(i,2)*B(2,j) + &
-! A2(i,3)*B(3,j) + &
-! A2(i,4)*B(4,j) + &
-! A2(i,5)*B(5,j)
-!
-! C3(i,j) = A3(i,1)*B(1,j) + &
-! A3(i,2)*B(2,j) + &
-! A3(i,3)*B(3,j) + &
-! A3(i,4)*B(4,j) + &
-! A3(i,5)*B(5,j)
-!
-! enddo
-! enddo
-!
-! end subroutine old_mxm_m2_m1_5points
Modified: seismo/3D/SPECFEM3D/trunk/compute_forces_elastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_forces_elastic_noDev.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/compute_forces_elastic_noDev.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -34,23 +34,18 @@
one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy,&
- epsilondev_xz,epsilondev_yz,iflag_attenuation_store,&
- rho_vs,&
+ epsilondev_xz,epsilondev_yz,epsilon_trace_over_3, &
+ iflag_attenuation_store,rho_vs,&
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,NGLOB_ADJOINT,NSPEC_ADJOINT, &
- b_displ,b_accel,kappa_kl,mu_kl,deltat, &
- NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ATT_AND_KERNEL,&
+ SIMULATION_TYPE,COMPUTE_AND_STORE_STRAIN,NSPEC_STRAIN_ONLY, &
+ NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ADJOINT,&
is_moho_top,is_moho_bot, &
- dsdx_top,dsdx_bot,b_dsdx_top,b_dsdx_bot, &
+ dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
- b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
- b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
- b_epsilondev_xz,b_epsilondev_yz, &
- b_alphaval,b_betaval,b_gammaval,&
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic,&
phase_ispec_inner_elastic )
@@ -88,14 +83,18 @@
! memory variables and standard linear solids for attenuation
logical :: ATTENUATION,USE_OLSEN_ATTENUATION
+ logical :: COMPUTE_AND_STORE_STRAIN
+ integer :: NSPEC_STRAIN_ONLY, NSPEC_ADJOINT
integer :: NSPEC_ATTENUATION_AB
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: iflag_attenuation_store
real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION) :: one_minus_sum_beta
real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: factor_common, alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS) :: &
R_xx,R_yy,R_xy,R_xz,R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB) :: &
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: epsilon_trace_over_3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vs
@@ -116,38 +115,18 @@
! adjoint simulations
integer :: SIMULATION_TYPE
- integer :: NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ATT_AND_KERNEL
- integer :: NGLOB_ADJOINT,NSPEC_ADJOINT
+ integer :: NSPEC_BOUN,NSPEC2D_MOHO
! moho kernel
real(kind=CUSTOM_REAL),dimension(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO):: &
- dsdx_top,dsdx_bot,b_dsdx_top,b_dsdx_bot
+ dsdx_top,dsdx_bot
logical,dimension(NSPEC_BOUN) :: is_moho_top,is_moho_bot
integer :: ispec2D_moho_top, ispec2D_moho_bot
- ! adjoint memory variables
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS) :: &
- b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL) :: &
- b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz
- real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: b_alphaval,b_betaval,b_gammaval
-
- ! adjoint wavefields
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT):: b_displ,b_accel
- ! adjoint kernels
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: &
- mu_kl, kappa_kl
- real(kind=CUSTOM_REAL) :: deltat
-
-!adjoint
-
! local parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
- integer ispec,iglob,ispec_p,num_elements
- integer i,j,k,l
-
real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
@@ -175,31 +154,13 @@
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
real(kind=CUSTOM_REAL) factor_loc,alphaval_loc,betaval_loc,gammaval_loc,Sn,Snp1
- real(kind=CUSTOM_REAL) epsilon_trace_over_3
+ real(kind=CUSTOM_REAL) templ
real(kind=CUSTOM_REAL) vs_val
integer i_SLS,iselected
+ integer ispec,iglob,ispec_p,num_elements
+ integer i,j,k,l
- ! adjoint backward arrays
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
- b_tempx1,b_tempx2,b_tempx3,b_tempy1,b_tempy2,b_tempy3,b_tempz1,b_tempz2,b_tempz3
- real(kind=CUSTOM_REAL):: dsxx,dsxy,dsxz,dsyy,dsyz,dszz
- real(kind=CUSTOM_REAL):: b_duxdxl,b_duxdyl,b_duxdzl,b_duydxl,b_duydyl,b_duydzl,b_duzdxl,b_duzdyl,b_duzdzl
- real(kind=CUSTOM_REAL):: b_duxdxl_plus_duydyl,b_duxdxl_plus_duzdzl,b_duydyl_plus_duzdzl
- real(kind=CUSTOM_REAL):: b_duxdyl_plus_duydxl,b_duzdxl_plus_duxdzl,b_duzdyl_plus_duydzl
- real(kind=CUSTOM_REAL):: b_dsxx,b_dsxy,b_dsxz,b_dsyy,b_dsyz,b_dszz
- real(kind=CUSTOM_REAL):: b_sigma_xx,b_sigma_yy,b_sigma_zz,b_sigma_xy,b_sigma_xz,b_sigma_yz
- real(kind=CUSTOM_REAL):: kappa_k, mu_k
- real(kind=CUSTOM_REAL) b_tempx1l,b_tempx2l,b_tempx3l
- real(kind=CUSTOM_REAL) b_tempy1l,b_tempy2l,b_tempy3l
- real(kind=CUSTOM_REAL) b_tempz1l,b_tempz2l,b_tempz3l
- ! local adjoint attenuation
- real(kind=CUSTOM_REAL) b_alphaval_loc,b_betaval_loc,b_gammaval_loc,b_Sn,b_Snp1
- real(kind=CUSTOM_REAL) b_epsilon_trace_over_3
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: b_epsilondev_xx_loc, &
- b_epsilondev_yy_loc, b_epsilondev_xy_loc, b_epsilondev_xz_loc, b_epsilondev_yz_loc
- real(kind=CUSTOM_REAL) b_R_xx_val,b_R_yy_val
- ! adjoint
if( iphase == 1 ) then
num_elements = nspec_outer_elastic
@@ -209,646 +170,345 @@
do ispec_p = 1,num_elements
-! if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+ ispec = phase_ispec_inner_elastic(ispec_p,iphase)
-! if( ispec_is_elastic(ispec) ) then
+ ! adjoint simulations: moho kernel
+ ! note: call this only once
+ if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
+ if (is_moho_top(ispec)) then
+ ispec2D_moho_top = ispec2D_moho_top + 1
+ else if (is_moho_bot(ispec)) then
+ ispec2D_moho_bot = ispec2D_moho_bot + 1
+ endif
+ endif
- ispec = phase_ispec_inner_elastic(ispec_p,iphase)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ tempx1l = 0.
+ tempx2l = 0.
+ tempx3l = 0.
- ! adjoint simulations: moho kernel
- if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
- if (is_moho_top(ispec)) then
- ispec2D_moho_top = ispec2D_moho_top + 1
- else if (is_moho_bot(ispec)) then
- ispec2D_moho_bot = ispec2D_moho_bot + 1
- endif
- endif
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ tempy1l = 0.
+ tempy2l = 0.
+ tempy3l = 0.
- tempx1l = 0.
- tempx2l = 0.
- tempx3l = 0.
+ tempz1l = 0.
+ tempz2l = 0.
+ tempz3l = 0.
- tempy1l = 0.
- tempy2l = 0.
- tempy3l = 0.
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempy1l = tempy1l + displ(2,iglob)*hp1
+ tempz1l = tempz1l + displ(3,iglob)*hp1
- tempz1l = 0.
- tempz2l = 0.
- tempz3l = 0.
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempy2l = tempy2l + displ(2,iglob)*hp2
+ tempz2l = tempz2l + displ(3,iglob)*hp2
- if (SIMULATION_TYPE == 3) then
- b_tempx1l = 0.
- b_tempx2l = 0.
- b_tempx3l = 0.
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3l = tempx3l + displ(1,iglob)*hp3
+ tempy3l = tempy3l + displ(2,iglob)*hp3
+ tempz3l = tempz3l + displ(3,iglob)*hp3
+ enddo
- b_tempy1l = 0.
- b_tempy2l = 0.
- b_tempy3l = 0.
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
- b_tempz1l = 0.
- b_tempz2l = 0.
- b_tempz3l = 0.
- endif
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
- do l=1,NGLLX
- hp1 = hprime_xx(i,l)
- iglob = ibool(l,j,k,ispec)
- tempx1l = tempx1l + displ(1,iglob)*hp1
- tempy1l = tempy1l + displ(2,iglob)*hp1
- tempz1l = tempz1l + displ(3,iglob)*hp1
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx1l = b_tempx1l + b_displ(1,iglob)*hp1
- b_tempy1l = b_tempy1l + b_displ(2,iglob)*hp1
- b_tempz1l = b_tempz1l + b_displ(3,iglob)*hp1
- endif ! adjoint
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- hp2 = hprime_yy(j,l)
- iglob = ibool(i,l,k,ispec)
- tempx2l = tempx2l + displ(1,iglob)*hp2
- tempy2l = tempy2l + displ(2,iglob)*hp2
- tempz2l = tempz2l + displ(3,iglob)*hp2
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx2l = b_tempx2l + b_displ(1,iglob)*hp2
- b_tempy2l = b_tempy2l + b_displ(2,iglob)*hp2
- b_tempz2l = b_tempz2l + b_displ(3,iglob)*hp2
- endif ! adjoint
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- hp3 = hprime_zz(k,l)
- iglob = ibool(i,j,l,ispec)
- tempx3l = tempx3l + displ(1,iglob)*hp3
- tempy3l = tempy3l + displ(2,iglob)*hp3
- tempz3l = tempz3l + displ(3,iglob)*hp3
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx3l = b_tempx3l + b_displ(1,iglob)*hp3
- b_tempy3l = b_tempy3l + b_displ(2,iglob)*hp3
- b_tempz3l = b_tempz3l + b_displ(3,iglob)*hp3
- endif ! adjoint
-
-
- enddo
+ ! adjoint simulations: save strain on the Moho boundary
+ if (SAVE_MOHO_MESH ) then
+ if (is_moho_top(ispec)) then
+ dsdx_top(1,1,i,j,k,ispec2D_moho_top) = duxdxl
+ dsdx_top(1,2,i,j,k,ispec2D_moho_top) = duxdyl
+ dsdx_top(1,3,i,j,k,ispec2D_moho_top) = duxdzl
+ dsdx_top(2,1,i,j,k,ispec2D_moho_top) = duydxl
+ dsdx_top(2,2,i,j,k,ispec2D_moho_top) = duydyl
+ dsdx_top(2,3,i,j,k,ispec2D_moho_top) = duydzl
+ dsdx_top(3,1,i,j,k,ispec2D_moho_top) = duzdxl
+ dsdx_top(3,2,i,j,k,ispec2D_moho_top) = duzdyl
+ dsdx_top(3,3,i,j,k,ispec2D_moho_top) = duzdzl
+ else if (is_moho_bot(ispec)) then
+ dsdx_bot(1,1,i,j,k,ispec2D_moho_bot) = duxdxl
+ dsdx_bot(1,2,i,j,k,ispec2D_moho_bot) = duxdyl
+ dsdx_bot(1,3,i,j,k,ispec2D_moho_bot) = duxdzl
+ dsdx_bot(2,1,i,j,k,ispec2D_moho_bot) = duydxl
+ dsdx_bot(2,2,i,j,k,ispec2D_moho_bot) = duydyl
+ dsdx_bot(2,3,i,j,k,ispec2D_moho_bot) = duydzl
+ dsdx_bot(3,1,i,j,k,ispec2D_moho_bot) = duzdxl
+ dsdx_bot(3,2,i,j,k,ispec2D_moho_bot) = duzdyl
+ dsdx_bot(3,3,i,j,k,ispec2D_moho_bot) = duzdzl
+ endif
+ endif
- ! get derivatives of ux, uy and uz with respect to x, y and z
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
- duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
- duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
- duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
- duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+ kappal = kappastore(i,j,k,ispec)
+ mul = mustore(i,j,k,ispec)
- duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
- duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
- duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+ if(ATTENUATION) then
+ ! uses scaling rule similar to Olsen et al. (2003) or mesh flag
+ if(USE_OLSEN_ATTENUATION) then
+ vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
+ call get_attenuation_model_olsen( vs_val, iselected )
+ else
+ ! iflag from (CUBIT) mesh
+ iselected = iflag_attenuation_store(i,j,k,ispec)
+ endif
- ! adjoint simulations: save strain on the Moho boundary
- if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
- if (is_moho_top(ispec)) then
- dsdx_top(1,1,i,j,k,ispec2D_moho_top) = duxdxl
- dsdx_top(1,2,i,j,k,ispec2D_moho_top) = duxdyl
- dsdx_top(1,3,i,j,k,ispec2D_moho_top) = duxdzl
- dsdx_top(2,1,i,j,k,ispec2D_moho_top) = duydxl
- dsdx_top(2,2,i,j,k,ispec2D_moho_top) = duydyl
- dsdx_top(2,3,i,j,k,ispec2D_moho_top) = duydzl
- dsdx_top(3,1,i,j,k,ispec2D_moho_top) = duzdxl
- dsdx_top(3,2,i,j,k,ispec2D_moho_top) = duzdyl
- dsdx_top(3,3,i,j,k,ispec2D_moho_top) = duzdzl
- else if (is_moho_bot(ispec)) then
- dsdx_bot(1,1,i,j,k,ispec2D_moho_bot) = duxdxl
- dsdx_bot(1,2,i,j,k,ispec2D_moho_bot) = duxdyl
- dsdx_bot(1,3,i,j,k,ispec2D_moho_bot) = duxdzl
- dsdx_bot(2,1,i,j,k,ispec2D_moho_bot) = duydxl
- dsdx_bot(2,2,i,j,k,ispec2D_moho_bot) = duydyl
- dsdx_bot(2,3,i,j,k,ispec2D_moho_bot) = duydzl
- dsdx_bot(3,1,i,j,k,ispec2D_moho_bot) = duzdxl
- dsdx_bot(3,2,i,j,k,ispec2D_moho_bot) = duzdyl
- dsdx_bot(3,3,i,j,k,ispec2D_moho_bot) = duzdzl
- endif
- endif
+ ! use unrelaxed parameters if attenuation
+ mul = mul * one_minus_sum_beta(iselected)
+
+ endif
- ! precompute some sums to save CPU time
- duxdxl_plus_duydyl = duxdxl + duydyl
- duxdxl_plus_duzdzl = duxdxl + duzdzl
- duydyl_plus_duzdzl = duydyl + duzdzl
- duxdyl_plus_duydxl = duxdyl + duydxl
- duzdxl_plus_duxdzl = duzdxl + duxdzl
- duzdyl_plus_duydzl = duzdyl + duydzl
+ ! full anisotropic case, stress calculations
+ if(ANISOTROPY) then
+ c11 = c11store(i,j,k,ispec)
+ c12 = c12store(i,j,k,ispec)
+ c13 = c13store(i,j,k,ispec)
+ c14 = c14store(i,j,k,ispec)
+ c15 = c15store(i,j,k,ispec)
+ c16 = c16store(i,j,k,ispec)
+ c22 = c22store(i,j,k,ispec)
+ c23 = c23store(i,j,k,ispec)
+ c24 = c24store(i,j,k,ispec)
+ c25 = c25store(i,j,k,ispec)
+ c26 = c26store(i,j,k,ispec)
+ c33 = c33store(i,j,k,ispec)
+ c34 = c34store(i,j,k,ispec)
+ c35 = c35store(i,j,k,ispec)
+ c36 = c36store(i,j,k,ispec)
+ c44 = c44store(i,j,k,ispec)
+ c45 = c45store(i,j,k,ispec)
+ c46 = c46store(i,j,k,ispec)
+ c55 = c55store(i,j,k,ispec)
+ c56 = c56store(i,j,k,ispec)
+ c66 = c66store(i,j,k,ispec)
- kappal = kappastore(i,j,k,ispec)
- mul = mustore(i,j,k,ispec)
+ sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+ sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+ sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+ sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+ sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+ sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- dsxx = duxdxl
- dsxy = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
- dsxz = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
- dsyy = duydyl
- dsyz = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
- dszz = duzdzl
+ else
- b_duxdxl = xixl*b_tempx1l + etaxl*b_tempx2l + gammaxl*b_tempx3l
- b_duxdyl = xiyl*b_tempx1l + etayl*b_tempx2l + gammayl*b_tempx3l
- b_duxdzl = xizl*b_tempx1l + etazl*b_tempx2l + gammazl*b_tempx3l
+ ! isotropic case
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
- b_duydxl = xixl*b_tempy1l + etaxl*b_tempy2l + gammaxl*b_tempy3l
- b_duydyl = xiyl*b_tempy1l + etayl*b_tempy2l + gammayl*b_tempy3l
- b_duydzl = xizl*b_tempy1l + etazl*b_tempy2l + gammazl*b_tempy3l
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
- b_duzdxl = xixl*b_tempz1l + etaxl*b_tempz2l + gammaxl*b_tempz3l
- b_duzdyl = xiyl*b_tempz1l + etayl*b_tempz2l + gammayl*b_tempz3l
- b_duzdzl = xizl*b_tempz1l + etazl*b_tempz2l + gammazl*b_tempz3l
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
- b_duxdxl_plus_duydyl = b_duxdxl + b_duydyl
- b_duxdxl_plus_duzdzl = b_duxdxl + b_duzdzl
- b_duydyl_plus_duzdzl = b_duydyl + b_duzdzl
- b_duxdyl_plus_duydxl = b_duxdyl + b_duydxl
- b_duzdxl_plus_duxdzl = b_duzdxl + b_duxdzl
- b_duzdyl_plus_duydzl = b_duzdyl + b_duydzl
+ endif ! ANISOTROPY
- b_dsxx = b_duxdxl
- b_dsxy = 0.5_CUSTOM_REAL * b_duxdyl_plus_duydxl
- b_dsxz = 0.5_CUSTOM_REAL * b_duzdxl_plus_duxdzl
- b_dsyy = b_duydyl
- b_dsyz = 0.5_CUSTOM_REAL * b_duzdyl_plus_duydzl
- b_dszz = b_duzdzl
+ ! subtract memory variables if attenuation
+ if(ATTENUATION) then
+ do i_sls = 1,N_SLS
+ R_xx_val = R_xx(i,j,k,ispec,i_sls)
+ R_yy_val = R_yy(i,j,k,ispec,i_sls)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
+ sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
+ sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
+ enddo
+ endif
- ! isotropic adjoint kernels: bulk (kappa) and shear (mu) kernels
- kappa_k = (duxdxl + duydyl + duzdzl) * (b_duxdxl + b_duydyl + b_duzdzl)
- mu_k = dsxx * b_dsxx + dsyy * b_dsyy + dszz * b_dszz + &
- 2._CUSTOM_REAL * (dsxy * b_dsxy + dsxz * b_dsxz + dsyz * b_dsyz) &
- - ONE_THIRD * kappa_k
-
- kappa_kl(i,j,k,ispec) = kappa_kl(i,j,k,ispec) + deltat * kappa_k
- mu_kl(i,j,k,ispec) = mu_kl(i,j,k,ispec) + 2._CUSTOM_REAL * deltat * mu_k
+ ! form dot product with test vector, symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
- if (SAVE_MOHO_MESH) then
- if (is_moho_top(ispec)) then
- b_dsdx_top(1,1,i,j,k,ispec2D_moho_top) = b_duxdxl
- b_dsdx_top(1,2,i,j,k,ispec2D_moho_top) = b_duxdyl
- b_dsdx_top(1,3,i,j,k,ispec2D_moho_top) = b_duxdzl
- b_dsdx_top(2,1,i,j,k,ispec2D_moho_top) = b_duydxl
- b_dsdx_top(2,2,i,j,k,ispec2D_moho_top) = b_duydyl
- b_dsdx_top(2,3,i,j,k,ispec2D_moho_top) = b_duydzl
- b_dsdx_top(3,1,i,j,k,ispec2D_moho_top) = b_duzdxl
- b_dsdx_top(3,2,i,j,k,ispec2D_moho_top) = b_duzdyl
- b_dsdx_top(3,3,i,j,k,ispec2D_moho_top) = b_duzdzl
- else if (is_moho_bot(ispec)) then
- b_dsdx_bot(1,1,i,j,k,ispec2D_moho_bot) = b_duxdxl
- b_dsdx_bot(1,2,i,j,k,ispec2D_moho_bot) = b_duxdyl
- b_dsdx_bot(1,3,i,j,k,ispec2D_moho_bot) = b_duxdzl
- b_dsdx_bot(2,1,i,j,k,ispec2D_moho_bot) = b_duydxl
- b_dsdx_bot(2,2,i,j,k,ispec2D_moho_bot) = b_duydyl
- b_dsdx_bot(2,3,i,j,k,ispec2D_moho_bot) = b_duydzl
- b_dsdx_bot(3,1,i,j,k,ispec2D_moho_bot) = b_duzdxl
- b_dsdx_bot(3,2,i,j,k,ispec2D_moho_bot) = b_duzdyl
- b_dsdx_bot(3,3,i,j,k,ispec2D_moho_bot) = b_duzdzl
- endif
- endif
- endif ! adjoint
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
- if(ATTENUATION) then
- ! compute deviatoric strain
- epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_xx_loc(i,j,k) = duxdxl - epsilon_trace_over_3
- epsilondev_yy_loc(i,j,k) = duydyl - epsilon_trace_over_3
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
-
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_epsilon_trace_over_3 = ONE_THIRD * (b_duxdxl + b_duydyl + b_duzdzl)
- b_epsilondev_xx_loc(i,j,k) = b_duxdxl - b_epsilon_trace_over_3
- b_epsilondev_yy_loc(i,j,k) = b_duydyl - b_epsilon_trace_over_3
- b_epsilondev_xy_loc(i,j,k) = 0.5 * b_duxdyl_plus_duydxl
- b_epsilondev_xz_loc(i,j,k) = 0.5 * b_duzdxl_plus_duxdzl
- b_epsilondev_yz_loc(i,j,k) = 0.5 * b_duzdyl_plus_duydzl
- endif ! adjoint
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
- ! uses scaling rule similar to Olsen et al. (2003) or mesh flag
- if(USE_OLSEN_ATTENUATION) then
- vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
- call get_attenuation_model_olsen( vs_val, iselected )
- else
- ! iflag from (CUBIT) mesh
- iselected = iflag_attenuation_store(i,j,k,ispec)
- endif
+ enddo
+ enddo
+ enddo
- ! use unrelaxed parameters if attenuation
- mul = mul * one_minus_sum_beta(iselected)
-
- endif
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
- ! full anisotropic case, stress calculations
- if(ANISOTROPY) then
- c11 = c11store(i,j,k,ispec)
- c12 = c12store(i,j,k,ispec)
- c13 = c13store(i,j,k,ispec)
- c14 = c14store(i,j,k,ispec)
- c15 = c15store(i,j,k,ispec)
- c16 = c16store(i,j,k,ispec)
- c22 = c22store(i,j,k,ispec)
- c23 = c23store(i,j,k,ispec)
- c24 = c24store(i,j,k,ispec)
- c25 = c25store(i,j,k,ispec)
- c26 = c26store(i,j,k,ispec)
- c33 = c33store(i,j,k,ispec)
- c34 = c34store(i,j,k,ispec)
- c35 = c35store(i,j,k,ispec)
- c36 = c36store(i,j,k,ispec)
- c44 = c44store(i,j,k,ispec)
- c45 = c45store(i,j,k,ispec)
- c46 = c46store(i,j,k,ispec)
- c55 = c55store(i,j,k,ispec)
- c56 = c56store(i,j,k,ispec)
- c66 = c66store(i,j,k,ispec)
- !if(ATTENUATION .and. not_fully_in_bedrock(ispec)) then
- ! mul = c44
- ! c11 = c11 + FOUR_THIRDS * minus_sum_beta * mul
- ! c12 = c12 - TWO_THIRDS * minus_sum_beta * mul
- ! c13 = c13 - TWO_THIRDS * minus_sum_beta * mul
- ! c22 = c22 + FOUR_THIRDS * minus_sum_beta * mul
- ! c23 = c23 - TWO_THIRDS * minus_sum_beta * mul
- ! c33 = c33 + FOUR_THIRDS * minus_sum_beta * mul
- ! c44 = c44 + minus_sum_beta * mul
- ! c55 = c55 + minus_sum_beta * mul
- ! c66 = c66 + minus_sum_beta * mul
- !endif
+ tempx1l = 0.
+ tempy1l = 0.
+ tempz1l = 0.
- sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
- c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
- sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
- c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
- sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
- c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
- sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
- c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
- sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
- c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
- sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
- c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
+ tempx2l = 0.
+ tempy2l = 0.
+ tempz2l = 0.
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_sigma_xx = c11*b_duxdxl + c16*b_duxdyl_plus_duydxl + c12*b_duydyl + &
- c15*b_duzdxl_plus_duxdzl + c14*b_duzdyl_plus_duydzl + c13*b_duzdzl
- b_sigma_yy = c12*b_duxdxl + c26*b_duxdyl_plus_duydxl + c22*b_duydyl + &
- c25*b_duzdxl_plus_duxdzl + c24*b_duzdyl_plus_duydzl + c23*b_duzdzl
- b_sigma_zz = c13*b_duxdxl + c36*b_duxdyl_plus_duydxl + c23*b_duydyl + &
- c35*b_duzdxl_plus_duxdzl + c34*b_duzdyl_plus_duydzl + c33*b_duzdzl
- b_sigma_xy = c16*b_duxdxl + c66*b_duxdyl_plus_duydxl + c26*b_duydyl + &
- c56*b_duzdxl_plus_duxdzl + c46*b_duzdyl_plus_duydzl + c36*b_duzdzl
- b_sigma_xz = c15*b_duxdxl + c56*b_duxdyl_plus_duydxl + c25*b_duydyl + &
- c55*b_duzdxl_plus_duxdzl + c45*b_duzdyl_plus_duydzl + c35*b_duzdzl
- b_sigma_yz = c14*b_duxdxl + c46*b_duxdyl_plus_duydxl + c24*b_duydyl + &
- c45*b_duzdxl_plus_duxdzl + c44*b_duzdyl_plus_duydzl + c34*b_duzdzl
- endif ! adjoint
- else
+ tempx3l = 0.
+ tempy3l = 0.
+ tempz3l = 0.
- ! isotropic case
- lambdalplus2mul = kappal + FOUR_THIRDS * mul
- lambdal = lambdalplus2mul - 2.*mul
+ do l=1,NGLLX
+ fac1 = hprimewgll_xx(l,i)
+ tempx1l = tempx1l + tempx1(l,j,k)*fac1
+ tempy1l = tempy1l + tempy1(l,j,k)*fac1
+ tempz1l = tempz1l + tempz1(l,j,k)*fac1
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
- sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
- sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ fac2 = hprimewgll_yy(l,j)
+ tempx2l = tempx2l + tempx2(i,l,k)*fac2
+ tempy2l = tempy2l + tempy2(i,l,k)*fac2
+ tempz2l = tempz2l + tempz2(i,l,k)*fac2
- sigma_xy = mul*duxdyl_plus_duydxl
- sigma_xz = mul*duzdxl_plus_duxdzl
- sigma_yz = mul*duzdyl_plus_duydzl
-
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_sigma_xx = lambdalplus2mul*b_duxdxl + lambdal*b_duydyl_plus_duzdzl
- b_sigma_yy = lambdalplus2mul*b_duydyl + lambdal*b_duxdxl_plus_duzdzl
- b_sigma_zz = lambdalplus2mul*b_duzdzl + lambdal*b_duxdxl_plus_duydyl
- b_sigma_xy = mul*b_duxdyl_plus_duydxl
- b_sigma_xz = mul*b_duzdxl_plus_duxdzl
- b_sigma_yz = mul*b_duzdyl_plus_duydzl
- endif !adjoint
-
- endif ! ANISOTROPY
-
- ! subtract memory variables if attenuation
- if(ATTENUATION) then
- do i_sls = 1,N_SLS
- R_xx_val = R_xx(i,j,k,ispec,i_sls)
- R_yy_val = R_yy(i,j,k,ispec,i_sls)
- sigma_xx = sigma_xx - R_xx_val
- sigma_yy = sigma_yy - R_yy_val
- sigma_zz = sigma_zz + R_xx_val + R_yy_val
- sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
- sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
- sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
-
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_R_xx_val = b_R_xx(i,j,k,ispec,i_sls)
- b_R_yy_val = b_R_yy(i,j,k,ispec,i_sls)
- b_sigma_xx = b_sigma_xx - b_R_xx_val
- b_sigma_yy = b_sigma_yy - b_R_yy_val
- b_sigma_zz = b_sigma_zz + b_R_xx_val + b_R_yy_val
- b_sigma_xy = b_sigma_xy - b_R_xy(i,j,k,ispec,i_sls)
- b_sigma_xz = b_sigma_xz - b_R_xz(i,j,k,ispec,i_sls)
- b_sigma_yz = b_sigma_yz - b_R_yz(i,j,k,ispec,i_sls)
- endif !adjoint
- enddo
- endif
-
- ! form dot product with test vector, symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
-
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx1(i,j,k) = jacobianl * (b_sigma_xx*xixl + b_sigma_xy*xiyl + b_sigma_xz*xizl)
- b_tempy1(i,j,k) = jacobianl * (b_sigma_xy*xixl + b_sigma_yy*xiyl + b_sigma_yz*xizl)
- b_tempz1(i,j,k) = jacobianl * (b_sigma_xz*xixl + b_sigma_yz*xiyl + b_sigma_zz*xizl)
- b_tempx2(i,j,k) = jacobianl * (b_sigma_xx*etaxl + b_sigma_xy*etayl + b_sigma_xz*etazl)
- b_tempy2(i,j,k) = jacobianl * (b_sigma_xy*etaxl + b_sigma_yy*etayl + b_sigma_yz*etazl)
- b_tempz2(i,j,k) = jacobianl * (b_sigma_xz*etaxl + b_sigma_yz*etayl + b_sigma_zz*etazl)
- b_tempx3(i,j,k) = jacobianl * (b_sigma_xx*gammaxl + b_sigma_xy*gammayl + b_sigma_xz*gammazl)
- b_tempy3(i,j,k) = jacobianl * (b_sigma_xy*gammaxl + b_sigma_yy*gammayl + b_sigma_yz*gammazl)
- b_tempz3(i,j,k) = jacobianl * (b_sigma_xz*gammaxl + b_sigma_yz*gammayl + b_sigma_zz*gammazl)
- endif !adjoint
-
- enddo
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ fac3 = hprimewgll_zz(l,k)
+ tempx3l = tempx3l + tempx3(i,j,l)*fac3
+ tempy3l = tempy3l + tempy3(i,j,l)*fac3
+ tempz3l = tempz3l + tempz3(i,j,l)*fac3
enddo
- enddo
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
- tempx1l = 0.
- tempy1l = 0.
- tempz1l = 0.
+ ! sum contributions from each element to the global mesh
+ iglob = ibool(i,j,k,ispec)
- tempx2l = 0.
- tempy2l = 0.
- tempz2l = 0.
+ accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
+ accel(2,iglob) = accel(2,iglob) - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
+ accel(3,iglob) = accel(3,iglob) - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
- tempx3l = 0.
- tempy3l = 0.
- tempz3l = 0.
+ ! update memory variables based upon the Runge-Kutta scheme
+ if(ATTENUATION) then
+
+ ! use Runge-Kutta scheme to march in time
+ do i_sls = 1,N_SLS
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx1l = 0.
- b_tempy1l = 0.
- b_tempz1l = 0.
- b_tempx2l = 0.
- b_tempy2l = 0.
- b_tempz2l = 0.
- b_tempx3l = 0.
- b_tempy3l = 0.
- b_tempz3l = 0.
- endif !adjoint
+ ! get coefficients for that standard linear solid
+ if( USE_OLSEN_ATTENUATION ) then
+ vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
+ call get_attenuation_model_olsen( vs_val, iselected )
+ else
+ iselected = iflag_attenuation_store(i,j,k,ispec)
+ endif
+
+ factor_loc = mustore(i,j,k,ispec) * factor_common(iselected,i_sls)
+
+ alphaval_loc = alphaval(iselected,i_sls)
+ betaval_loc = betaval(iselected,i_sls)
+ gammaval_loc = gammaval(iselected,i_sls)
+
+ ! term in xx
+ Sn = factor_loc * epsilondev_xx(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xx_loc(i,j,k)
+ R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
- do l=1,NGLLX
- fac1 = hprimewgll_xx(l,i)
- tempx1l = tempx1l + tempx1(l,j,k)*fac1
- tempy1l = tempy1l + tempy1(l,j,k)*fac1
- tempz1l = tempz1l + tempz1(l,j,k)*fac1
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx1l = b_tempx1l + b_tempx1(l,j,k)*fac1
- b_tempy1l = b_tempy1l + b_tempy1(l,j,k)*fac1
- b_tempz1l = b_tempz1l + b_tempz1(l,j,k)*fac1
- endif
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+ ! term in yy
+ Sn = factor_loc * epsilondev_yy(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_yy_loc(i,j,k)
+ R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- fac2 = hprimewgll_yy(l,j)
- tempx2l = tempx2l + tempx2(i,l,k)*fac2
- tempy2l = tempy2l + tempy2(i,l,k)*fac2
- tempz2l = tempz2l + tempz2(i,l,k)*fac2
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx2l = b_tempx2l + b_tempx2(i,l,k)*fac2
- b_tempy2l = b_tempy2l + b_tempy2(i,l,k)*fac2
- b_tempz2l = b_tempz2l + b_tempz2(i,l,k)*fac2
- endif
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+ ! term in zz not computed since zero trace
+
+ ! term in xy
+ Sn = factor_loc * epsilondev_xy(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xy_loc(i,j,k)
+ R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
+
+ ! term in xz
+ Sn = factor_loc * epsilondev_xz(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xz_loc(i,j,k)
+ R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- fac3 = hprimewgll_zz(l,k)
- tempx3l = tempx3l + tempx3(i,j,l)*fac3
- tempy3l = tempy3l + tempy3(i,j,l)*fac3
- tempz3l = tempz3l + tempz3(i,j,l)*fac3
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_tempx3l = b_tempx3l + b_tempx3(i,j,l)*fac3
- b_tempy3l = b_tempy3l + b_tempy3(i,j,l)*fac3
- b_tempz3l = b_tempz3l + b_tempz3(i,j,l)*fac3
- endif
- enddo
+ ! term in yz
+ Sn = factor_loc * epsilondev_yz(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_yz_loc(i,j,k)
+ R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
+
+ enddo ! end of loop on memory variables
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
+ endif ! end attenuation
- ! sum contributions from each element to the global mesh
- iglob = ibool(i,j,k,ispec)
-
- accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
- accel(2,iglob) = accel(2,iglob) - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
- accel(3,iglob) = accel(3,iglob) - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
-
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_accel(1,iglob) = b_accel(1,iglob) - (fac1*b_tempx1l + fac2*b_tempx2l + fac3*b_tempx3l)
- b_accel(2,iglob) = b_accel(2,iglob) - (fac1*b_tempy1l + fac2*b_tempy2l + fac3*b_tempy3l)
- b_accel(3,iglob) = b_accel(3,iglob) - (fac1*b_tempz1l + fac2*b_tempz2l + fac3*b_tempz3l)
- endif !adjoint
-
- ! update memory variables based upon the Runge-Kutta scheme
- if(ATTENUATION) then
-
- ! use Runge-Kutta scheme to march in time
- do i_sls = 1,N_SLS
-
- ! get coefficients for that standard linear solid
- if( USE_OLSEN_ATTENUATION ) then
- vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
- call get_attenuation_model_olsen( vs_val, iselected )
- else
- iselected = iflag_attenuation_store(i,j,k,ispec)
- endif
-
- factor_loc = mustore(i,j,k,ispec) * factor_common(iselected,i_sls)
-
- alphaval_loc = alphaval(iselected,i_sls)
- betaval_loc = betaval(iselected,i_sls)
- gammaval_loc = gammaval(iselected,i_sls)
-
- ! term in xx
- Sn = factor_loc * epsilondev_xx(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xx_loc(i,j,k)
- R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
-
- ! term in yy
- Sn = factor_loc * epsilondev_yy(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_yy_loc(i,j,k)
- R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
-
- ! term in zz not computed since zero trace
-
- ! term in xy
- Sn = factor_loc * epsilondev_xy(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xy_loc(i,j,k)
- R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
-
- ! term in xz
- Sn = factor_loc * epsilondev_xz(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xz_loc(i,j,k)
- R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
-
- ! term in yz
- Sn = factor_loc * epsilondev_yz(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_yz_loc(i,j,k)
- R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
-
- !adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_alphaval_loc = b_alphaval(iselected,i_sls)
- b_betaval_loc = b_betaval(iselected,i_sls)
- b_gammaval_loc = b_gammaval(iselected,i_sls)
- ! term in xx
- b_Sn = factor_loc * b_epsilondev_xx(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_xx_loc(i,j,k)
- b_R_xx(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xx(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- ! term in yy
- b_Sn = factor_loc * b_epsilondev_yy(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_yy_loc(i,j,k)
- b_R_yy(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_yy(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- ! term in zz not computed since zero trace
- ! term in xy
- b_Sn = factor_loc * b_epsilondev_xy(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_xy_loc(i,j,k)
- b_R_xy(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xy(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- ! term in xz
- b_Sn = factor_loc * b_epsilondev_xz(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_xz_loc(i,j,k)
- b_R_xz(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xz(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- ! term in yz
- b_Sn = factor_loc * b_epsilondev_yz(i,j,k,ispec)
- b_Snp1 = factor_loc * b_epsilondev_yz_loc(i,j,k)
- b_R_yz(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_yz(i,j,k,ispec,i_sls) + &
- b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
- endif !adjoint
-
- enddo ! end of loop on memory variables
-
- endif ! end attenuation
-
-
- enddo
- enddo
enddo
+ enddo
+ enddo
- ! save deviatoric strain for Runge-Kutta scheme
- if(ATTENUATION) then
- epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
- epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
- epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
- epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
- epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) then
- b_epsilondev_xx(:,:,:,ispec) = b_epsilondev_xx_loc(:,:,:)
- b_epsilondev_yy(:,:,:,ispec) = b_epsilondev_yy_loc(:,:,:)
- b_epsilondev_xy(:,:,:,ispec) = b_epsilondev_xy_loc(:,:,:)
- b_epsilondev_xz(:,:,:,ispec) = b_epsilondev_xz_loc(:,:,:)
- b_epsilondev_yz(:,:,:,ispec) = b_epsilondev_yz_loc(:,:,:)
- endif !adjoint
- endif
-! endif ! ispec_is_elastic
-! endif ! if (ispec_is_inner(ispec) .eqv. phase_is_inner)
+ ! save deviatoric strain for Runge-Kutta scheme
+ if ( COMPUTE_AND_STORE_STRAIN ) then
+ epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
+ epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
+ epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
+ epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
+ epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
+ endif
enddo ! spectral element loop
-! forces in elastic media calculated in compute_forces_elastic...
-!! adding source
-! do isource = 1,NSOURCES
-!
-! if (ispec_is_inner(ispec_selected_source(isource)) .eqv. phase_is_inner) then
-!
-! if(USE_FORCE_POINT_SOURCE) then
-!
-!! add the source (only if this proc carries the source)
-! if(myrank == islice_selected_source(isource)) then
-!
-! iglob = ibool(nint(xi_source(isource)), &
-! nint(eta_source(isource)), &
-! nint(gamma_source(isource)), &
-! ispec_selected_source(isource))
-! f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
-! t0 = 1.2d0/f0
-!
-! if (it == 1 .and. myrank == 0) then
-! print *,'using a source of dominant frequency ',f0
-! print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
-! print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
-! endif
-!
-! ! we use nu_source(:,3) here because we want a source normal to the surface.
-! ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
-! !accel(:,iglob) = accel(:,iglob) + &
-! ! sngl(nu_source(:,3,isource) * 10000000.d0 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
-! ! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
-! accel(:,iglob) = accel(:,iglob) + &
-! sngl(nu_source(:,3,isource) * 1.d10 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
-! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
-!
-! endif
-! endif
-!
-! endif
-!
-! enddo
end subroutine compute_forces_elastic_noDev
Modified: seismo/3D/SPECFEM3D/trunk/compute_interpolated_dva.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_interpolated_dva.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/compute_interpolated_dva.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -125,7 +125,7 @@
dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
! acoustic elements
-! returns displacement/velocity/acceleration (dxd,..,vxd,..,axd,.. ) at receiver location
+! returns displacement/velocity/pressure (dxd,..,vxd,..,axd,.. ) at receiver location
implicit none
include 'constants.h'
Added: seismo/3D/SPECFEM3D/trunk/compute_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_kernels.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/compute_kernels.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -0,0 +1,141 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine compute_kernels()
+
+! kernel calculations
+! see e.g. Tromp et al. (2005)
+
+ use specfem_par
+ use specfem_par_elastic
+ use specfem_par_acoustic
+
+ implicit none
+ ! local parameters
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ):: b_displ_elm,accel_elm
+ real(kind=CUSTOM_REAL) :: kappal
+ integer :: i,j,k,ispec,iglob
+ real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc,b_epsilondev_loc
+
+ ! updates kernels
+ do ispec = 1, NSPEC_AB
+
+ ! elastic domains
+ if( ispec_is_elastic(ispec) ) then
+
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+
+ ! isotropic kernels
+ ! note: takes displacement from backward/reconstructed (forward) field b_displ
+ ! and acceleration from adjoint field accel (containing adjoint sources)
+ !
+ ! note: : time integral summation uses deltat
+ !
+ ! compare with Tromp et al. (2005), eq. (14), which takes adjoint displacement
+ ! and forward acceleration, that is the symmetric form of what is calculated here
+ ! however, this kernel expression is symmetric with regards
+ ! to interchange adjoint - forward field
+ rho_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) &
+ + deltat * dot_product(accel(:,iglob), b_displ(:,iglob))
+
+ ! kernel for shear modulus, see e.g. Tromp et al. (2005), equation (17)
+ ! note: multiplication with 2*mu(x) will be done after the time loop
+ epsilondev_loc(1) = epsilondev_xx(i,j,k,ispec)
+ epsilondev_loc(2) = epsilondev_yy(i,j,k,ispec)
+ epsilondev_loc(3) = epsilondev_xy(i,j,k,ispec)
+ epsilondev_loc(4) = epsilondev_xz(i,j,k,ispec)
+ epsilondev_loc(5) = epsilondev_yz(i,j,k,ispec)
+
+ b_epsilondev_loc(1) = b_epsilondev_xx(i,j,k,ispec)
+ b_epsilondev_loc(2) = b_epsilondev_yy(i,j,k,ispec)
+ b_epsilondev_loc(3) = b_epsilondev_xy(i,j,k,ispec)
+ b_epsilondev_loc(4) = b_epsilondev_xz(i,j,k,ispec)
+ b_epsilondev_loc(5) = b_epsilondev_yz(i,j,k,ispec)
+
+ mu_kl(i,j,k,ispec) = mu_kl(i,j,k,ispec) &
+ + deltat * (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
+ + (epsilondev_loc(1)+epsilondev_loc(2)) * (b_epsilondev_loc(1)+b_epsilondev_loc(2)) &
+ + 2 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) + &
+ epsilondev_loc(5)*b_epsilondev_loc(5)) )
+
+ ! kernel for bulk modulus, see e.g. Tromp et al. (2005), equation (18)
+ ! note: multiplication with kappa(x) will be done after the time loop
+ kappa_kl(i,j,k,ispec) = kappa_kl(i,j,k,ispec) &
+ + deltat * (9 * epsilon_trace_over_3(i,j,k,ispec) &
+ * b_epsilon_trace_over_3(i,j,k,ispec))
+
+ enddo
+ enddo
+ enddo
+ endif !ispec_is_elastic
+
+ ! acoustic domains
+ if( ispec_is_acoustic(ispec) ) then
+
+ ! backward fields: displacement vector
+ call compute_gradient(ispec,NSPEC_ADJOINT,NGLOB_ADJOINT, &
+ b_potential_acoustic, b_displ_elm,&
+ hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ ibool,rhostore)
+ ! adjoint fields: acceleration vector
+ call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
+ potential_dot_dot_acoustic, accel_elm,&
+ hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ ibool,rhostore)
+
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+
+ ! density kernel
+ rho_ac_kl(i,j,k,ispec) = rho_ac_kl(i,j,k,ispec) &
+ - deltat * dot_product(accel_elm(:,i,j,k), b_displ_elm(:,i,j,k))
+
+ ! bulk modulus kernel
+ kappal = kappastore(i,j,k,ispec)
+ kappa_ac_kl(i,j,k,ispec) = kappa_ac_kl(i,j,k,ispec) &
+ - deltat / kappal &
+ * potential_dot_dot_acoustic(iglob) &
+ * b_potential_dot_dot_acoustic(iglob)
+
+ enddo
+ enddo
+ enddo
+ endif ! ispec_is_acoustic
+
+ enddo
+
+ ! moho kernel
+ if( ELASTIC_SIMULATION .and. SAVE_MOHO_MESH ) then
+ call compute_boundary_kernel()
+ endif
+
+ end subroutine compute_kernels
Modified: seismo/3D/SPECFEM3D/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/constants.h.in 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/constants.h.in 2010-08-29 15:46:15 UTC (rev 17141)
@@ -121,12 +121,12 @@
! number of steps to save the state variables in the forward simulation,
! to be used in the backward reconstruction in the presence of attenuation
- integer, parameter :: NSTEP_Q_SAVE = 200
+ integer, parameter :: NSTEP_Q_SAVE = 50 ! depending on stability of reconstruction, up to 200
! the scratch disk to save the state variables saved in the forward
! simulation, this can be a global scratch disk in case you run out of
! space on the local scratch disk
- character(len=256), parameter :: LOCAL_PATH_Q = '/ibrixfs1/scratch/lqy/DATABASES_MPI_Q/'
+ character(len=256), parameter :: LOCAL_PATH_Q = 'DATABASES_MPI' ! or e.g. '/ibrixfs1/scratch/lqy/DATABASES_MPI_Q/'
!------------------------------------------------------
! nlegoff -- Variables that should be read/computed elsewhere.
@@ -140,7 +140,7 @@
! this can be useful e.g. for oil industry foothills simulations or asteroid simulations
! in which the source is a vertical force, normal force, impact etc.
logical, parameter :: USE_FORCE_POINT_SOURCE = .false.
- double precision, parameter :: FACTOR_FORCE_SOURCE = 1.d10
+ double precision, parameter :: FACTOR_FORCE_SOURCE = 1.d15
! the receivers can be located inside the model
logical, parameter :: RECVS_CAN_BE_BURIED_EXT_MESH = .true.
Modified: seismo/3D/SPECFEM3D/trunk/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/finalize_simulation.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/finalize_simulation.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -51,20 +51,21 @@
write(27) displ
write(27) veloc
write(27) accel
+
+ if (ATTENUATION) then
+ write(27) R_xx
+ write(27) R_yy
+ write(27) R_xy
+ write(27) R_xz
+ write(27) R_yz
+ write(27) epsilondev_xx
+ write(27) epsilondev_yy
+ write(27) epsilondev_xy
+ write(27) epsilondev_xz
+ write(27) epsilondev_yz
+ endif
endif
- if (ATTENUATION) then
- write(27) R_xx
- write(27) R_yy
- write(27) R_xy
- write(27) R_xz
- write(27) R_yz
- write(27) epsilondev_xx
- write(27) epsilondev_yy
- write(27) epsilondev_xy
- write(27) epsilondev_xz
- write(27) epsilondev_yz
- endif
close(27)
! adjoint simulations
Modified: seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -114,6 +114,14 @@
! if attenuation is off, set dummy size of arrays to one
NSPEC_ATTENUATION_AB = 1
endif
+ ! needed for attenuation and/or kernel computations
+ if( ATTENUATION .or. SIMULATION_TYPE == 3 ) then
+ COMPUTE_AND_STORE_STRAIN = .true.
+ NSPEC_STRAIN_ONLY = NSPEC_AB
+ else
+ COMPUTE_AND_STORE_STRAIN = .false.
+ NSPEC_STRAIN_ONLY = 1
+ endif
! anisotropy arrays size
if( ANISOTROPY ) then
Modified: seismo/3D/SPECFEM3D/trunk/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/iterate_time.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/iterate_time.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -94,7 +94,7 @@
! adjoint simulations: kernels
if( SIMULATION_TYPE == 3 ) then
- call it_update_adjointkernels()
+ call compute_kernels()
endif
! outputs movie files
@@ -354,6 +354,12 @@
if( it > 1 .and. it < NSTEP) then
! adjoint simulations
+
+! note backward/reconstructed wavefields:
+! storing wavefield displ() at time step it, corresponds to time (it-1)*DT - t0 (see routine write_seismograms_to_file )
+! reconstucted wavefield b_displ() at it corresponds to time (NSTEP-it-1)*DT - t0
+! we read in the reconstructed wavefield at the end of the time iteration loop, i.e. after the Newark scheme,
+! thus, indexing is NSTEP-it (rather than something like NSTEP-(it-1) )
if (SIMULATION_TYPE == 3 .and. mod(NSTEP-it,NSTEP_Q_SAVE) == 0) then
! reads files content
write(outputname,"('save_Q_arrays_',i6.6,'.bin')") NSTEP-it
@@ -401,109 +407,12 @@
write(27) epsilondev_yz
endif
if( ACOUSTIC_SIMULATION ) then
- write(27) b_potential_acoustic
- write(27) b_potential_dot_acoustic
- write(27) b_potential_dot_dot_acoustic
+ write(27) potential_acoustic
+ write(27) potential_dot_acoustic
+ write(27) potential_dot_dot_acoustic
endif
close(27)
endif ! SIMULATION_TYPE
endif ! it
end subroutine it_store_attenuation_arrays
-
-!================================================================
-
- subroutine it_update_adjointkernels()
-
-! kernel calculations
-
- use specfem_par
- use specfem_par_elastic
- use specfem_par_acoustic
-
- implicit none
- ! local parameters
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ):: b_displ_elm,accel_elm
- real(kind=CUSTOM_REAL) :: kappal
- integer :: i,j,k,ispec,iglob
-
- !elastic domains
- if(ELASTIC_SIMULATION ) then
-
- ! NOTE: kappa and mu kernels have already been updated in compute_forces_elastic()
-
- ! density kernel update
- do ispec = 1, NSPEC_AB
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
-
- ! note: takes displacement from backward/reconstructed (forward) field b_displ
- ! and acceleration from adjoint field accel (containing adjoint sources)
- !
- ! note: : time integral summation uses deltat
- !
- ! compare with Tromp et al. (2005), eq. (14), which takes adjoint displacement
- ! and forward acceleration, that is the symmetric form of what is calculated here
- ! however, this kernel expression is symmetric with regards to interchange adjoint - forward field
- rho_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) &
- + deltat * dot_product(accel(:,iglob), b_displ(:,iglob))
- enddo
- enddo
- enddo
- enddo
-
- ! moho kernel
- if (SAVE_MOHO_MESH) then
- call compute_boundary_kernel()
- endif
-
- endif ! elastic
-
- ! acoustic domains
- if( ACOUSTIC_SIMULATION ) then
-
- do ispec=1,NSPEC_AB
-
- ! acoustic wave field
- if( ispec_is_acoustic(ispec) ) then
-
- ! backward fields: displacement vector
- call compute_gradient(ispec,NSPEC_ADJOINT,NGLOB_ADJOINT, &
- b_potential_acoustic, b_displ_elm,&
- hprime_xx,hprime_yy,hprime_zz, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- ibool,rhostore)
- ! adjoint fields: acceleration vector
- call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
- potential_dot_dot_acoustic, accel_elm,&
- hprime_xx,hprime_yy,hprime_zz, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- ibool,rhostore)
-
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
-
- ! density kernel
- rho_ac_kl(i,j,k,ispec) = rho_ac_kl(i,j,k,ispec) &
- - deltat * dot_product(accel_elm(:,i,j,k), b_displ_elm(:,i,j,k))
-
- ! bulk modulus kernel
- kappal = kappastore(i,j,k,ispec)
- kappa_ac_kl(i,j,k,ispec) = kappa_ac_kl(i,j,k,ispec) &
- - deltat * kappal &
- * potential_dot_dot_acoustic(iglob)/kappal &
- * b_potential_dot_dot_acoustic(iglob)/kappal
- enddo
- enddo
- enddo
-
- endif ! ispec_is_acoustic
- enddo
- endif !acoustic
-
- end subroutine it_update_adjointkernels
-
Modified: seismo/3D/SPECFEM3D/trunk/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/prepare_timerun.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/prepare_timerun.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -348,26 +348,8 @@
+ deltat**4*tauinv(:,:)**3 / 24._CUSTOM_REAL
gammaval(:,:) = deltat / 2._CUSTOM_REAL + deltat**2*tauinv(:,:) / 6._CUSTOM_REAL &
+ deltat**3*tauinv(:,:)**2 / 24._CUSTOM_REAL
- endif
-
- !pll, to put elsewhere
- ! note: currently, they need to be defined here, as they are used in the routine arguments
- ! for compute_forces_with_Deville()
- allocate(R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(epsilondev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
- allocate(epsilondev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
- allocate(epsilondev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
- allocate(epsilondev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
- allocate(epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
-
-! clear memory variables if attenuation
- if(ATTENUATION) then
-
+ ! clear memory variables if attenuation
! initialize memory variables for attenuation
epsilondev_xx(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_yy(:,:,:,:) = 0._CUSTOM_REAL
@@ -478,22 +460,23 @@
read(27) b_displ
read(27) b_veloc
read(27) b_accel
- endif
-
- ! memory variables if attenuation
- if( ATTENUATION ) then
- read(27) b_R_xx
- read(27) b_R_yy
- read(27) b_R_xy
- read(27) b_R_xz
- read(27) b_R_yz
- read(27) b_epsilondev_xx
- read(27) b_epsilondev_yy
- read(27) b_epsilondev_xy
- read(27) b_epsilondev_xz
- read(27) b_epsilondev_yz
- endif
+ ! memory variables if attenuation
+ if( ATTENUATION ) then
+ read(27) b_R_xx
+ read(27) b_R_yy
+ read(27) b_R_xy
+ read(27) b_R_xz
+ read(27) b_R_yz
+ read(27) b_epsilondev_xx
+ read(27) b_epsilondev_yy
+ read(27) b_epsilondev_xy
+ read(27) b_epsilondev_xz
+ read(27) b_epsilondev_yz
+ endif
+
+ endif
+
close(27)
endif
Modified: seismo/3D/SPECFEM3D/trunk/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/read_mesh_databases.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/read_mesh_databases.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -128,6 +128,24 @@
allocate(c56store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
allocate(c66store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ ! note: currently, they need to be defined, as they are used in the routine arguments
+ ! for compute_forces_elastic_Deville()
+ allocate(R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+
+ ! needed for attenuation and/or kernel computations
+ allocate(epsilondev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
+ allocate(epsilondev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
+ allocate(epsilondev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
+ allocate(epsilondev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
+ allocate(epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
+
+ ! note: needed for argument of deville routine
+ allocate(epsilon_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+
read(27) rmass
if( OCEANS ) then
! ocean mass matrix
@@ -177,6 +195,7 @@
call exit_MPI(myrank,'error attenuation flag entry exceeds range')
endif
endif
+
! absorbing boundary surface
read(27) num_abs_boundary_faces
@@ -432,6 +451,24 @@
allocate(b_request_recv_vector_ext_mesh(num_interfaces_ext_mesh))
allocate(b_buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
allocate(b_buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
+
+ ! allocates attenuation solids
+ if( ATTENUATION ) then
+ allocate(b_R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
+ b_R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
+ b_R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
+ b_R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
+ b_R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS) )
+ endif
+
+ ! note: these arrays are needed for attenuation and/or kernel computations
+ allocate(b_epsilondev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
+ b_epsilondev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
+ b_epsilondev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
+ b_epsilondev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
+ b_epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY) )
+ ! needed for kernel computations
+ allocate(b_epsilon_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
endif
@@ -455,21 +492,6 @@
allocate(b_buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
endif
-
-! allocates attenuation solids
- if( ATTENUATION .and. SIMULATION_TYPE == 3 ) then
- allocate(b_R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS) )
-
- allocate(b_epsilondev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL), &
- b_epsilondev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL), &
- b_epsilondev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL), &
- b_epsilondev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL), &
- b_epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL) )
- endif
! ADJOINT moho
! moho boundary
Modified: seismo/3D/SPECFEM3D/trunk/save_adjoint_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/save_adjoint_kernels.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/save_adjoint_kernels.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -33,6 +33,7 @@
implicit none
integer:: ispec,i,j,k,iglob
+ real(kind=CUSTOM_REAL) :: rhol,mul,kappal
! finalizes calculation of rhop, beta, alpha kernels
do ispec = 1, NSPEC_AB
@@ -46,29 +47,30 @@
iglob = ibool(i,j,k,ispec)
! isotropic adjoint kernels (see e.g. Tromp et al. 2005)
+ rhol = rho_vs(i,j,k,ispec)**2 / mustore(i,j,k,ispec)
+ mul = mustore(i,j,k,ispec)
+ kappal = kappastore(i,j,k,ispec)
! density kernel
! multiplies with rho
- rho_kl(i,j,k,ispec) = - rho_vs(i,j,k,ispec)**2 / mustore(i,j,k,ispec) * rho_kl(i,j,k,ispec)
+ rho_kl(i,j,k,ispec) = - rhol * rho_kl(i,j,k,ispec)
! shear modulus kernel
- mu_kl(i,j,k,ispec) = - mustore(i,j,k,ispec) * mu_kl(i,j,k,ispec)
+ mu_kl(i,j,k,ispec) = - mul * mu_kl(i,j,k,ispec)
! bulk modulus kernel
- kappa_kl(i,j,k,ispec) = - kappastore(i,j,k,ispec) * kappa_kl(i,j,k,ispec)
+ kappa_kl(i,j,k,ispec) = - kappal * kappa_kl(i,j,k,ispec)
! density prime kernel
rhop_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) + kappa_kl(i,j,k,ispec) + mu_kl(i,j,k,ispec)
! vs kernel
beta_kl(i,j,k,ispec) = 2._CUSTOM_REAL * (mu_kl(i,j,k,ispec) &
- - 4._CUSTOM_REAL * mustore(i,j,k,ispec) &
- / (3._CUSTOM_REAL * kappastore(i,j,k,ispec)) * kappa_kl(i,j,k,ispec))
+ - 4._CUSTOM_REAL * mul / (3._CUSTOM_REAL * kappal) * kappa_kl(i,j,k,ispec))
! vp kernel
alpha_kl(i,j,k,ispec) = 2._CUSTOM_REAL * (1._CUSTOM_REAL &
- + 4._CUSTOM_REAL * mustore(i,j,k,ispec) &
- / (3._CUSTOM_REAL * kappastore(i,j,k,ispec))) * kappa_kl(i,j,k,ispec)
+ + 4._CUSTOM_REAL * mul / (3._CUSTOM_REAL * kappal) ) * kappa_kl(i,j,k,ispec)
enddo
enddo
enddo
@@ -84,7 +86,7 @@
! rho prime kernel
rhop_ac_kl(i,j,k,ispec) = rho_ac_kl(i,j,k,ispec) + kappa_ac_kl(i,j,k,ispec)
- ! vp kernel
+ ! kappa kernel
alpha_ac_kl(i,j,k,ispec) = TWO * kappa_ac_kl(i,j,k,ispec)
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -124,7 +124,9 @@
hdur_gaussian = hdur/SOURCE_DECAY_MIMIC_TRIANGLE
! define t0 as the earliest start time
- t0 = - 1.5d0 * minval(t_cmt-hdur)
+ ! note: an earlier start time also reduces numerical noise due to a
+ ! non-zero offset at the beginning of the source time function
+ t0 = - 2.0d0 * minval(t_cmt-hdur) ! - 1.5d0 * minval(t_cmt-hdur)
! uses an earlier start time if source is acoustic with a gaussian source time function
t0_ac = 0.0d0
Modified: seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -259,6 +259,7 @@
R_xx,R_yy,R_xy,R_xz,R_yz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: epsilon_trace_over_3
! displacement, velocity, acceleration
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: displ,veloc,accel
@@ -277,6 +278,10 @@
c55store,c56store,c66store
integer :: NSPEC_ANISO
+! for attenuation and/or kernel simulations
+ integer :: NSPEC_STRAIN_ONLY
+ logical :: COMPUTE_AND_STORE_STRAIN
+
! material flag
logical, dimension(:), allocatable :: ispec_is_elastic
integer, dimension(:,:), allocatable :: phase_ispec_inner_elastic
@@ -297,6 +302,8 @@
b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_epsilon_trace_over_3
+
integer:: NSPEC_ATT_AND_KERNEL
! adjoint kernels
Modified: seismo/3D/SPECFEM3D/trunk/write_PNM_GIF_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/write_PNM_GIF_data.f90 2010-08-28 22:47:36 UTC (rev 17140)
+++ seismo/3D/SPECFEM3D/trunk/write_PNM_GIF_data.f90 2010-08-29 15:46:15 UTC (rev 17141)
@@ -200,7 +200,7 @@
if( count /= num_iglob_image_surface) call exit_mpi(myrank,'error image point number')
- !daniel: outputs found global points into vtk file
+ !daniel: outputs global points into vtk file
!vtkfilename = prname(1:len_trim(prname))//'GIF_image_points'
!call write_VTK_data_points(NGLOB_AB,xstore,ystore,zstore, &
! iglob_coord,count,vtkfilename)
More information about the CIG-COMMITS
mailing list