[cig-commits] r15666 - in seismo/2D/SPECFEM2D/trunk: . DATA
cmorency at geodynamics.org
cmorency at geodynamics.org
Mon Sep 14 14:56:42 PDT 2009
Author: cmorency
Date: 2009-09-14 14:56:41 -0700 (Mon, 14 Sep 2009)
New Revision: 15666
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt
seismo/2D/SPECFEM2D/trunk/adj_seismogram.f90
seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90
seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/plotpost.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
Added SH (membrane) waves calculation. Also fixed some bugs related to acoustic kernels.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2009-09-14 21:56:41 UTC (rev 15666)
@@ -30,6 +30,7 @@
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
Q0 = 1 # quality factor for viscous attenuation
freq0 = 10 # frequency for viscous attenuation
+body_waves = .true. # set the type of calculation (P-SV or SH/membrane waves)
# absorbing boundaries parameters
absorbing_conditions = .true. # absorbing boundary active or not
Modified: seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt
===================================================================
--- seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt 2009-09-14 21:56:41 UTC (rev 15666)
@@ -126,7 +126,7 @@
isolver/save_forward combination and a warning message appears in the ouput
file)
-Important output files (for example, for the elastic case)
+Important output files (for example, for the elastic case, P-SV waves)
absorb_elastic_bottom*****.bin
absorb_elastic_left*****.bin
absorb_elastic_right*****.bin
@@ -136,13 +136,15 @@
S****.AA.BHZ.semd
Second: define the adjoint source
-Use adj_seismogram.f90 which is in UTILS/adjoint.
+Use adj_seismogram.f90
Edit to update NSTEP, nrec, t0, deltat, and the position of the cut to pic
any given phase if needed (tstart,tend), add the right number of stations, and
put one component of the source to zero if needed.
-The ouput files of adj_seismogram.f90 are S****.AA.BHX.adj and S****.AA.BHZ.adj. They need to be
-kept in the OUTPUT_FILES directory together with the absorb_elastic_****.bin
-and lastframe_elastic.bin files to be read when running the adjoint simulation.
+The ouput files of adj_seismogram.f90 are S****.AA.BHX.adj and S****.AA.BHZ.adj, for P-SV waves (and
+S****.AA.BHY.adj, for SH (membrane) waves). Note that you will need these three
+files (S****.AA.BHX.adj, S****.AA.BHY.adj and S****.AA.BHZ.adj) to be present in the OUTPUT_FILES directory
+together with the absorb_elastic_****.bin and lastframe_elastic.bin files to be read
+when running the adjoint simulation.
Third: run the adjoint simulation
Make sure that the adjoint source files absorbing boundaries and last frame files are
@@ -155,7 +157,13 @@
snapshot_rhop_alpha_beta*****
which are the primary moduli kernels and the phase velocities kernels respectively.
-Note: At the moment, adjoint simulations do not support anisotropy, attenuation, and viscous damping.
+Note1: At the moment, adjoint simulations do not support anisotropy, attenuation, and viscous damping.
+Note2: You will need S****.AA.BHX.adj, S****.AA.BHY.adj and S****.AA.BHZ.adj
+to be present in OUTPUT_FILES even if you are just running an acoustic or
+poroelastic adjoint simulation.
+S****.AA.BHX.adj is the only relevant component for an acoustic case.
+S****.AA.BHX.adj and S****.AA.BHZ.adj are the only relevant components for a
+poroelastic case.
--------------------------------------------------
COUPLED SIMULATIONS
@@ -167,6 +175,19 @@
Elastic/poroelastic coupling supports anisotropy, but not attenuation for the
elastic material.
+
+How to run P-SV or SH (membrane) waves simulation :
+---------------------------------------------------
+To run a P-SV waves calculation propagating in the x-z plane,
+set body_waves = .true. in the Par_file.
+To run a SH (membrane) waves calculation traveling in the x-z plane with a
+y-component of motion, set body_waves = .false.
+
+This feature is only implemented for elastic materials and sensitivity kernels
+can be calculated (see Tape, Liu & Tromp, GJI 2006 for details on membrane
+surface waves).
+
+
--------------------------
--------------------------
--------------------------
Modified: seismo/2D/SPECFEM2D/trunk/adj_seismogram.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/adj_seismogram.f90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/adj_seismogram.f90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -5,37 +5,45 @@
implicit none
!
-! user edit
- integer, parameter :: NSTEP = 6000
+!!!! user edit
+ integer, parameter :: NSTEP = 26000
integer, parameter :: nrec = 1
- double precision, parameter :: t0 = 0.4
- double precision, parameter :: deltat = 1d-3
+ double precision, parameter :: t0 = 8
+ double precision, parameter :: deltat = 0.25d-2
double precision, parameter :: EPS = 1.d-40
-!
- integer :: itime,icomp,istart,iend,nlen,irec
+!!!!
+ integer :: itime,icomp,istart,iend,nlen,irec,NDIM,NDIMr,adj_comp
double precision :: time,tstart(nrec),tend(nrec)
character(len=150), dimension(nrec) :: station_name
double precision, dimension(NSTEP) :: time_window
- double precision :: seism(NSTEP,2),Nnorm,seism_win(NSTEP)
+ double precision :: seism(NSTEP,3),Nnorm,seism_win(NSTEP)
double precision :: seism_veloc(NSTEP),seism_accel(NSTEP),ft_bar(NSTEP)
- character(len=3) :: comp(2)
+ character(len=3) :: compr(2),comp(3)
character(len=150) :: filename,filename2
- include 'constants.h'
+ NDIM=3
+ comp = (/"BHX","BHY","BHZ"/)
-! user edit
+!!!! user edit
+! which calculation: P-SV (use (1)) or SH (membrane) (use (2)) waves
+! NDIMr=2 !(1)
+ NDIMr=1 !(2)
+! list of stations
station_name(1) = 'S0001'
- tstart(1) = 3.5d0 + t0
- tend(1) = 4.3d0 + t0
-!
+ tstart(1) = 39.7d0 + t0
+ tend(1) = 41d0 + t0
+! which calculation: P-SV (use (1)) or SH (membrane) (use (2)) waves
+! compr = (/"BHX","BHZ"/) !(1)
+ compr = (/"BHY","dummy"/) !(2)
+! chose the component for the adjoint source (adj_comp = 1: X, 2:Y, 3:Z)
+ adj_comp = 2
+!!!!
- comp = (/"BHX","BHZ"/)
-
do irec =1,nrec
- do icomp = 1, NDIM
+ do icomp = 1, NDIMr
- filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// comp(icomp) // '.semd'
+ filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// compr(icomp) // '.semd'
open(unit = 10, file = trim(filename))
do itime = 1,NSTEP
@@ -44,9 +52,18 @@
enddo
+ if(NDIMr==2)then
+ seism(:,3) = seism(:,2)
+ seism(:,2) = 0.d0
+ else
+ seism(:,2) = seism(:,1)
+ seism(:,1) = 0.d0
+ seism(:,3) = 0.d0
+ endif
+
close(10)
+
-
istart = max(floor(tstart(irec)/deltat),1)
iend = min(floor(tend(irec)/deltat),NSTEP)
print*,'istart =',istart, 'iend =', iend
@@ -56,6 +73,8 @@
do icomp = 1, NDIM
+ print*,comp(icomp)
+
filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// comp(icomp) // '.adj'
open(unit = 11, file = trim(filename))
@@ -85,8 +104,8 @@
! Nnorm = deltat * sum(time_window(:) * seism_veloc(:) * seism_veloc(:))
! cross-correlation traveltime adjoint source
if(abs(Nnorm) > EPS) then
- ft_bar(:) = - seism_veloc(:) * time_window(:) / Nnorm
-! ft_bar(:) = seism_veloc(:) * time_window(:) / Nnorm
+! ft_bar(:) = - seism_veloc(:) * time_window(:) / Nnorm
+ ft_bar(:) = seism_veloc(:) * time_window(:) / Nnorm
print*,'Norm =', Nnorm
else
print *, 'norm < EPS for file '
@@ -94,9 +113,8 @@
ft_bar(:) = 0.d0
endif
-! user edit: which component
do itime =1,NSTEP
- if(icomp == 1) then
+ if(icomp == adj_comp) then
write(11,*) (itime-1)*deltat - t0, ft_bar(itime)
else
write(11,*) (itime-1)*deltat - t0, 0.d0
@@ -107,6 +125,8 @@
close(11)
enddo
+ print*,'*************************'
+ print*,'The input files (S****.AA.BHX/BHY/BHZ.adj) needed to run the adjoint simulation are in OUTPUT_FILES'
+ print*,'*************************'
-
end program adj_seismogram
Modified: seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -145,7 +145,7 @@
character(len=*) adj_source_file
! output
- real(kind=CUSTOM_REAL), dimension(NSTEP,NDIM,NGLLX,NGLLZ) :: adj_sourcearray
+ real(kind=CUSTOM_REAL), dimension(NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearray
! Gauss-Lobatto-Legendre points of integration and weights
double precision, dimension(NGLLX) :: xigll
@@ -153,11 +153,11 @@
double precision :: hxir(NGLLX), hpxir(NGLLX), hgammar(NGLLZ), hpgammar(NGLLZ)
- real(kind=CUSTOM_REAL) :: adj_src_s(NSTEP,NDIM)
+ real(kind=CUSTOM_REAL) :: adj_src_s(NSTEP,3)
integer icomp, itime, i, k, ios
double precision :: junk
- character(len=3) :: comp(2)
+ character(len=3) :: comp(3)
character(len=150) :: filename
call lagrange_any(xi_receiver,NGLLX,xigll,hxir,hpxir)
@@ -165,9 +165,9 @@
adj_sourcearray(:,:,:,:) = 0.
- comp = (/"BHX","BHZ"/)
+ comp = (/"BHX","BHY","BHZ"/)
- do icomp = 1, NDIM
+ do icomp = 1,3
filename = 'OUTPUT_FILES/'//trim(adj_source_file) // '.'// comp(icomp) // '.adj'
open(unit = IIN, file = trim(filename), iostat = ios)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -53,7 +53,7 @@
jbegin_left,jend_left,jbegin_right,jend_right,isolver,save_forward,b_absorb_acoustic_left,&
b_absorb_acoustic_right,b_absorb_acoustic_bottom,&
b_absorb_acoustic_top,nspec_xmin,nspec_xmax,&
- nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax,kappa_ac_k)
+ nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax)
! compute forces for the acoustic elements
@@ -87,7 +87,6 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,rhoext
- real(kind=CUSTOM_REAL), dimension(npoin) :: kappa_ac_k
double precision, dimension(NGLLZ,nspec_xmin,NSTEP) :: b_absorb_acoustic_left
double precision, dimension(NGLLZ,nspec_xmax,NSTEP) :: b_absorb_acoustic_right
double precision, dimension(NGLLX,nspec_zmax,NSTEP) :: b_absorb_acoustic_top
@@ -170,15 +169,11 @@
dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
- if(isolver == 2) then
+ if(isolver == 2) then
b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
+ endif
-! kernels calculation
- iglob = ibool(i,j,ispec)
- kappa_ac_k(iglob) = dux_dxl * b_dux_dxl + dux_dzl * b_dux_dzl
- endif
-
jacobianl = jacobian(i,j,ispec)
! if external density model
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -42,7 +42,7 @@
!
!========================================================================
- subroutine compute_forces_elastic(npoin,nspec,myrank,nelemabs,numat, &
+ subroutine compute_forces_elastic(body_waves,npoin,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
@@ -66,6 +66,7 @@
include "constants.h"
+ logical :: body_waves
integer :: NSOURCE, i_source
integer :: npoin,nspec,myrank,nelemabs,numat,it,NSTEP
integer, dimension(NSOURCE) :: ispec_selected_source,is_proc_source,source_type
@@ -92,7 +93,8 @@
logical, dimension(nspec) :: elastic
logical, dimension(4,nelemabs) :: codeabs
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: accel_elastic,veloc_elastic,displ_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin) :: accel_elastic,veloc_elastic,displ_elastic
+ real(kind=CUSTOM_REAL), dimension(2,npoin) :: displ_att
double precision, dimension(2,numat) :: density
double precision, dimension(4,3,numat) :: elastcoef
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
@@ -100,13 +102,13 @@
real(kind=CUSTOM_REAL), dimension(NSOURCE,NSTEP) :: source_time_function
real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: b_accel_elastic,b_displ_elastic
- real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,NDIM,NGLLX,NGLLZ) :: adj_sourcearrays
+ real(kind=CUSTOM_REAL), dimension(3,npoin) :: b_accel_elastic,b_displ_elastic
+ real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(npoin) :: mu_k,kappa_k
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLZ,nspec_xmin,NSTEP) :: b_absorb_elastic_left
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLZ,nspec_xmax,NSTEP) :: b_absorb_elastic_right
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,nspec_zmax,NSTEP) :: b_absorb_elastic_top
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,nspec_zmin,NSTEP) :: b_absorb_elastic_bottom
+ real(kind=CUSTOM_REAL), dimension(3,NGLLZ,nspec_xmin,NSTEP) :: b_absorb_elastic_left
+ real(kind=CUSTOM_REAL), dimension(3,NGLLZ,nspec_xmax,NSTEP) :: b_absorb_elastic_right
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,nspec_zmax,NSTEP) :: b_absorb_elastic_top
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,nspec_zmin,NSTEP) :: b_absorb_elastic_bottom
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11,e13
@@ -134,18 +136,18 @@
integer :: ispec,i,j,k,iglob,ispecabs,ibegin,iend,irec,irec_local
! spatial derivatives
- real(kind=CUSTOM_REAL) :: dux_dxi,dux_dgamma,duz_dxi,duz_dgamma
- real(kind=CUSTOM_REAL) :: dux_dxl,duz_dxl,dux_dzl,duz_dzl
- real(kind=CUSTOM_REAL) :: b_dux_dxi,b_dux_dgamma,b_duz_dxi,b_duz_dgamma
- real(kind=CUSTOM_REAL) :: b_dux_dxl,b_duz_dxl,b_dux_dzl,b_duz_dzl
+ real(kind=CUSTOM_REAL) :: dux_dxi,dux_dgamma,duy_dxi,duy_dgamma,duz_dxi,duz_dgamma
+ real(kind=CUSTOM_REAL) :: dux_dxl,duy_dxl,duz_dxl,dux_dzl,duy_dzl,duz_dzl
+ real(kind=CUSTOM_REAL) :: b_dux_dxi,b_dux_dgamma,b_duy_dxi,b_duy_dgamma,b_duz_dxi,b_duz_dgamma
+ real(kind=CUSTOM_REAL) :: b_dux_dxl,b_duy_dxl,b_duz_dxl,b_dux_dzl,b_duy_dzl,b_duz_dzl
real(kind=CUSTOM_REAL) :: dsxx,dsxz,dszz
real(kind=CUSTOM_REAL) :: b_dsxx,b_dsxz,b_dszz
- real(kind=CUSTOM_REAL) :: sigma_xx,sigma_xz,sigma_zz
- real(kind=CUSTOM_REAL) :: b_sigma_xx,b_sigma_xz,b_sigma_zz
- real(kind=CUSTOM_REAL) :: nx,nz,vx,vz,vn,rho_vp,rho_vs,tx,tz,weight,xxi,zxi,xgamma,zgamma,jacobian1D
+ real(kind=CUSTOM_REAL) :: sigma_xx,sigma_xy,sigma_xz,sigma_zy,sigma_zz
+ real(kind=CUSTOM_REAL) :: b_sigma_xx,b_sigma_xy,b_sigma_xz,b_sigma_zy,b_sigma_zz
+ real(kind=CUSTOM_REAL) :: nx,nz,vx,vy,vz,vn,rho_vp,rho_vs,tx,ty,tz,weight,xxi,zxi,xgamma,zgamma,jacobian1D
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: tempx1,tempx2,tempz1,tempz2
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: b_tempx1,b_tempx2,b_tempz1,b_tempz2
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: tempx1,tempx2,tempy1,tempy2,tempz1,tempz2
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: b_tempx1,b_tempx2,b_tempy1,b_tempy2,b_tempz1,b_tempz2
! Jacobian matrix and determinant
real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl,jacobianl
@@ -173,8 +175,12 @@
integer :: ifirstelem,ilastelem
! compute Grad(displ_elastic) at time step n for attenuation
- if(TURN_ATTENUATION_ON) call compute_gradient_attenuation(displ_elastic,dux_dxl_n,duz_dxl_n, &
+ if(TURN_ATTENUATION_ON) then
+ displ_att(1,:) = displ_elastic(1,:)
+ displ_att(2,:) = displ_elastic(3,:)
+ call compute_gradient_attenuation(displ_att,dux_dxl_n,duz_dxl_n, &
dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,npoin)
+ endif
ifirstelem = 1
ilastelem = nspec
@@ -182,6 +188,21 @@
! loop over spectral elements
do ispec = ifirstelem,ilastelem
+ tempx1(:,:) = ZERO
+ tempy1(:,:) = ZERO
+ tempz1(:,:) = ZERO
+ tempx2(:,:) = ZERO
+ tempy2(:,:) = ZERO
+ tempz2(:,:) = ZERO
+ if(isolver ==2)then
+ b_tempx1(:,:) = ZERO
+ b_tempy1(:,:) = ZERO
+ b_tempz1(:,:) = ZERO
+ b_tempx2(:,:) = ZERO
+ b_tempy2(:,:) = ZERO
+ b_tempz2(:,:) = ZERO
+ endif
+
!---
!--- elastic spectral element
!---
@@ -208,16 +229,20 @@
! derivative along x and along z
dux_dxi = ZERO
+ duy_dxi = ZERO
duz_dxi = ZERO
dux_dgamma = ZERO
+ duy_dgamma = ZERO
duz_dgamma = ZERO
if(isolver == 2) then ! Adjoint calculation, backward wavefield
b_dux_dxi = ZERO
+ b_duy_dxi = ZERO
b_duz_dxi = ZERO
b_dux_dgamma = ZERO
+ b_duy_dgamma = ZERO
b_duz_dgamma = ZERO
endif
@@ -225,15 +250,19 @@
! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
- duz_dxi = duz_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ duy_dxi = duy_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ duz_dxi = duz_dxi + displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
- duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ duy_dgamma = duy_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
if(isolver == 2) then ! Adjoint calculation, backward wavefield
b_dux_dxi = b_dux_dxi + b_displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
- b_duz_dxi = b_duz_dxi + b_displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ b_duy_dxi = b_duy_dxi + b_displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ b_duz_dxi = b_duz_dxi + b_displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
b_dux_dgamma = b_dux_dgamma + b_displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
- b_duz_dgamma = b_duz_dgamma + b_displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ b_duy_dgamma = b_duy_dgamma + b_displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ b_duz_dgamma = b_duz_dgamma + b_displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
endif
enddo
@@ -246,6 +275,9 @@
dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+ duy_dxl = duy_dxi*xixl + duy_dgamma*gammaxl
+ duy_dzl = duy_dxi*xizl + duy_dgamma*gammazl
+
duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
@@ -253,6 +285,9 @@
b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
+ b_duy_dxl = b_duy_dxi*xixl + b_duy_dgamma*gammaxl
+ b_duy_dzl = b_duy_dxi*xizl + b_duy_dgamma*gammazl
+
b_duz_dxl = b_duz_dxi*xixl + b_duz_dgamma*gammaxl
b_duz_dzl = b_duz_dxi*xizl + b_duz_dgamma*gammazl
endif
@@ -297,12 +332,16 @@
! no attenuation
sigma_xx = lambdalplus2mul_relaxed*dux_dxl + lambdal_relaxed*duz_dzl
+ sigma_xy = mul_relaxed*duy_dxl
sigma_xz = mul_relaxed*(duz_dxl + dux_dzl)
+ sigma_zy = mul_relaxed*duy_dzl
sigma_zz = lambdalplus2mul_relaxed*duz_dzl + lambdal_relaxed*dux_dxl
if(isolver == 2) then ! Adjoint calculation, backward wavefield
b_sigma_xx = lambdalplus2mul_relaxed*b_dux_dxl + lambdal_relaxed*b_duz_dzl
+ b_sigma_xy = mul_relaxed*b_duy_dxl
b_sigma_xz = mul_relaxed*(b_duz_dxl + b_dux_dzl)
+ b_sigma_zy = mul_relaxed*b_duy_dzl
b_sigma_zz = lambdalplus2mul_relaxed*b_duz_dzl + lambdal_relaxed*b_dux_dxl
endif
@@ -321,6 +360,7 @@
! Pre-kernels calculation
if(isolver == 2) then
iglob = ibool(i,j,ispec)
+ if(body_waves)then !P-SV waves
dsxx = dux_dxl
dsxz = HALF * (duz_dxl + dux_dzl)
dszz = duz_dzl
@@ -332,6 +372,9 @@
kappa_k(iglob) = (dux_dxl + duz_dzl) * (b_dux_dxl + b_duz_dzl)
mu_k(iglob) = dsxx * b_dsxx + dszz * b_dszz + &
2._CUSTOM_REAL * dsxz * b_dsxz - 1._CUSTOM_REAL/3._CUSTOM_REAL * kappa_k(iglob)
+ else !SH (membrane) waves
+ mu_k(iglob) = duy_dxl * b_duy_dxl + duy_dzl * b_duy_dzl
+ endif
endif
jacobianl = jacobian(i,j,ispec)
@@ -339,16 +382,20 @@
! weak formulation term based on stress tensor (non-symmetric form)
! also add GLL integration weights
tempx1(i,j) = wzgll(j)*jacobianl*(sigma_xx*xixl+sigma_xz*xizl)
+ tempy1(i,j) = wzgll(j)*jacobianl*(sigma_xy*xixl+sigma_zy*xizl)
tempz1(i,j) = wzgll(j)*jacobianl*(sigma_xz*xixl+sigma_zz*xizl)
tempx2(i,j) = wxgll(i)*jacobianl*(sigma_xx*gammaxl+sigma_xz*gammazl)
+ tempy2(i,j) = wxgll(i)*jacobianl*(sigma_xy*gammaxl+sigma_zy*gammazl)
tempz2(i,j) = wxgll(i)*jacobianl*(sigma_xz*gammaxl+sigma_zz*gammazl)
if(isolver == 2) then ! Adjoint calculation, backward wavefield
b_tempx1(i,j) = wzgll(j)*jacobianl*(b_sigma_xx*xixl+b_sigma_xz*xizl)
+ b_tempy1(i,j) = wzgll(j)*jacobianl*(b_sigma_xy*xixl+b_sigma_zy*xizl)
b_tempz1(i,j) = wzgll(j)*jacobianl*(b_sigma_xz*xixl+b_sigma_zz*xizl)
b_tempx2(i,j) = wxgll(i)*jacobianl*(b_sigma_xx*gammaxl+b_sigma_xz*gammazl)
+ b_tempy2(i,j) = wxgll(i)*jacobianl*(b_sigma_xy*gammaxl+b_sigma_zy*gammazl)
b_tempz2(i,j) = wxgll(i)*jacobianl*(b_sigma_xz*gammaxl+b_sigma_zz*gammazl)
endif
@@ -368,12 +415,15 @@
! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tempx1(k,j)*hprimewgll_xx(k,i) + tempx2(i,k)*hprimewgll_zz(k,j))
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tempz1(k,j)*hprimewgll_xx(k,i) + tempz2(i,k)*hprimewgll_zz(k,j))
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tempy1(k,j)*hprimewgll_xx(k,i) + tempy2(i,k)*hprimewgll_zz(k,j))
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tempz1(k,j)*hprimewgll_xx(k,i) + tempz2(i,k)*hprimewgll_zz(k,j))
if(isolver == 2) then ! Adjoint calculation, backward wavefield
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
(b_tempx1(k,j)*hprimewgll_xx(k,i) + b_tempx2(i,k)*hprimewgll_zz(k,j))
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
+ (b_tempy1(k,j)*hprimewgll_xx(k,i) + b_tempy2(i,k)*hprimewgll_zz(k,j))
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
(b_tempz1(k,j)*hprimewgll_xx(k,i) + b_tempz2(i,k)*hprimewgll_zz(k,j))
endif
enddo
@@ -459,22 +509,33 @@
! Clayton-Engquist condition if elastic
if(elastic(ispec)) then
vx = veloc_elastic(1,iglob) - veloc_horiz
- vz = veloc_elastic(2,iglob) - veloc_vert
+ vy = veloc_elastic(2,iglob)
+ vz = veloc_elastic(3,iglob) - veloc_vert
vn = nx*vx+nz*vz
tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tz + traction_z_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
if(save_forward .and. isolver ==1) then
+ if(body_waves)then !P-SV waves
b_absorb_elastic_left(1,j,ib_xmin(ispecabs),it) = tx*weight
- b_absorb_elastic_left(2,j,ib_xmin(ispecabs),it) = tz*weight
+ b_absorb_elastic_left(3,j,ib_xmin(ispecabs),it) = tz*weight
+ else !SH (membrane) waves
+ b_absorb_elastic_left(2,j,ib_xmin(ispecabs),it) = ty*weight
+ endif
elseif(isolver == 2) then
+ if(body_waves)then !P-SV waves
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_left(1,j,ib_xmin(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_left(3,j,ib_xmin(ispecabs),NSTEP-it+1)
+ else !SH (membrane) waves
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_left(2,j,ib_xmin(ispecabs),NSTEP-it+1)
+ endif
endif
endif
@@ -536,23 +597,33 @@
! Clayton-Engquist condition if elastic
if(elastic(ispec)) then
vx = veloc_elastic(1,iglob) - veloc_horiz
- vz = veloc_elastic(2,iglob) - veloc_vert
+ vy = veloc_elastic(2,iglob)
+ vz = veloc_elastic(3,iglob) - veloc_vert
vn = nx*vx+nz*vz
tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tz - traction_z_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight
-
if(save_forward .and. isolver ==1) then
+ if(body_waves)then !P-SV waves
b_absorb_elastic_right(1,j,ib_xmax(ispecabs),it) = tx*weight
- b_absorb_elastic_right(2,j,ib_xmax(ispecabs),it) = tz*weight
+ b_absorb_elastic_right(3,j,ib_xmax(ispecabs),it) = tz*weight
+ else! SH (membrane) waves
+ b_absorb_elastic_right(2,j,ib_xmax(ispecabs),it) = ty*weight
+ endif
elseif(isolver == 2) then
+ if(body_waves)then !P-SV waves
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_right(1,j,ib_xmax(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_right(3,j,ib_xmax(ispecabs),NSTEP-it+1)
+ else! SH (membrane) waves
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_right(2,j,ib_xmax(ispecabs),NSTEP-it+1)
+ endif
endif
endif
@@ -620,22 +691,33 @@
! Clayton-Engquist condition if elastic
if(elastic(ispec)) then
vx = veloc_elastic(1,iglob) - veloc_horiz
- vz = veloc_elastic(2,iglob) - veloc_vert
+ vy = veloc_elastic(2,iglob)
+ vz = veloc_elastic(3,iglob) - veloc_vert
vn = nx*vx+nz*vz
tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tz + traction_z_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
if(save_forward .and. isolver ==1) then
+ if(body_waves)then !P-SV waves
b_absorb_elastic_bottom(1,i,ib_zmin(ispecabs),it) = tx*weight
- b_absorb_elastic_bottom(2,i,ib_zmin(ispecabs),it) = tz*weight
+ b_absorb_elastic_bottom(3,i,ib_zmin(ispecabs),it) = tz*weight
+ else!SH (membrane) waves
+ b_absorb_elastic_bottom(2,i,ib_zmin(ispecabs),it) = ty*weight
+ endif
elseif(isolver == 2) then
+ if(body_waves)then !P-SV waves
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_bottom(1,i,ib_zmin(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_bottom(3,i,ib_zmin(ispecabs),NSTEP-it+1)
+ else!SH (membrane) waves
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_bottom(2,i,ib_zmin(ispecabs),NSTEP-it+1)
+ endif
endif
endif
@@ -695,22 +777,33 @@
! Clayton-Engquist condition if elastic
if(elastic(ispec)) then
vx = veloc_elastic(1,iglob) - veloc_horiz
- vz = veloc_elastic(2,iglob) - veloc_vert
+ vy = veloc_elastic(2,iglob)
+ vz = veloc_elastic(3,iglob) - veloc_vert
vn = nx*vx+nz*vz
tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tz - traction_z_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight
if(save_forward .and. isolver ==1) then
+ if(body_waves)then !P-SV waves
b_absorb_elastic_top(1,i,ib_zmax(ispecabs),it) = tx*weight
- b_absorb_elastic_top(2,i,ib_zmax(ispecabs),it) = tz*weight
+ b_absorb_elastic_top(3,i,ib_zmax(ispecabs),it) = tz*weight
+ else!SH (membrane) waves
+ b_absorb_elastic_top(2,i,ib_zmax(ispecabs),it) = ty*weight
+ endif
elseif(isolver == 2) then
+ if(body_waves)then !P-SV waves
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_top(1,i,ib_zmax(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_top(3,i,ib_zmax(ispecabs),NSTEP-it+1)
+ else!SH (membrane) waves
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_top(2,i,ib_zmax(ispecabs),NSTEP-it+1)
+ endif
endif
endif
@@ -733,26 +826,32 @@
! moment tensor
if(source_type(i_source) == 2) then
+ if(.not.body_waves) call exit_MPI('cannot have moment tensor source in SH (membrane) waves calculation')
+
if(isolver == 1) then ! forward wavefield
! add source array
do j=1,NGLLZ
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_source(i_source))
- accel_elastic(:,iglob) = accel_elastic(:,iglob) + &
- sourcearray(i_source,:,i,j)*source_time_function(i_source,it)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + &
+ sourcearray(i_source,1,i,j)*source_time_function(i_source,it)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + &
+ sourcearray(i_source,2,i,j)*source_time_function(i_source,it)
enddo
enddo
else ! backward wavefield
do j=1,NGLLZ
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_source(i_source))
- b_accel_elastic(:,iglob) = b_accel_elastic(:,iglob) + &
- sourcearray(i_source,:,i,j)*source_time_function(i_source,NSTEP-it+1)
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + &
+ sourcearray(i_source,1,i,j)*source_time_function(i_source,NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) + &
+ sourcearray(i_source,2,i,j)*source_time_function(i_source,NSTEP-it+1)
enddo
enddo
endif !endif isolver == 1
- endif
+ endif !if(source_type(i_source) == 2)
endif ! if this processor carries the source and the source element is elastic
enddo ! do i_source=1,NSOURCE
@@ -769,7 +868,12 @@
do j=1,NGLLZ
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_rec(irec))
- accel_elastic(:,iglob) = accel_elastic(:,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,:,i,j)
+ if(body_waves)then !P-SH waves
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,3,i,j)
+ else !SH (membrane) waves
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,2,i,j)
+ endif
enddo
enddo
@@ -784,7 +888,7 @@
if(TURN_ATTENUATION_ON) then
! compute Grad(displ_elastic) at time step n+1 for attenuation
- call compute_gradient_attenuation(displ_elastic,dux_dxl_np1,duz_dxl_np1, &
+ call compute_gradient_attenuation(displ_att,dux_dxl_np1,duz_dxl_np1, &
dux_dzl_np1,duz_dzl_np1,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,npoin)
! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -103,7 +103,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
real(kind=CUSTOM_REAL), dimension(NSOURCE,NSTEP) :: source_time_function
real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
- real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,NDIM,NGLLX,NGLLZ) :: adj_sourcearrays
+ real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(npoin) :: C_k,M_k
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLZ,nspec_xmin,NSTEP) :: b_absorb_poro_w_left
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLZ,nspec_xmax,NSTEP) :: b_absorb_poro_w_right
@@ -847,8 +847,10 @@
do j=1,NGLLZ
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_rec(irec))
- accelw_poroelastic(:,iglob) = accelw_poroelastic(:,iglob) - &
- rhol_f/rhol_bar*adj_sourcearrays(irec_local,NSTEP-it+1,:,i,j)
+ accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) - &
+ rhol_f/rhol_bar*adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
+ accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) - &
+ rhol_f/rhol_bar*adj_sourcearrays(irec_local,NSTEP-it+1,3,i,j)
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -103,7 +103,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
real(kind=CUSTOM_REAL), dimension(NSOURCE,NSTEP) :: source_time_function
real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
- real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,NDIM,NGLLX,NGLLZ) :: adj_sourcearrays
+ real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(npoin) :: mufr_k,B_k
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLZ,nspec_xmin,NSTEP) :: b_absorb_poro_s_left
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLZ,nspec_xmax,NSTEP) :: b_absorb_poro_s_right
@@ -864,7 +864,8 @@
do j=1,NGLLZ
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_rec(irec))
- accels_poroelastic(:,iglob) = accels_poroelastic(:,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,:,i,j)
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,3,i,j)
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -68,8 +68,9 @@
logical, dimension(nspec) :: elastic,poroelastic
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_dot_acoustic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displ_elastic,displs_poroelastic,displw_poroelastic
- double precision, dimension(NDIM,npoin) :: vector_field_display
+ real(kind=CUSTOM_REAL), dimension(3,npoin) :: displ_elastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,displw_poroelastic
+ double precision, dimension(3,npoin) :: vector_field_display
! array with derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
@@ -101,7 +102,7 @@
do j = 1,NGLLZ
do i = 1,NGLLX
iglob = ibool(i,j,ispec)
- vector_field_display(2,iglob) = pressure_element(i,j)
+ vector_field_display(3,iglob) = pressure_element(i,j)
enddo
enddo
@@ -142,7 +143,8 @@
logical, dimension(nspec) :: elastic,poroelastic
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_dot_acoustic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displ_elastic,displs_poroelastic,displw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin) :: displ_elastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,displw_poroelastic
! array with derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
@@ -231,9 +233,9 @@
! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
- duz_dxi = duz_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ duz_dxi = duz_dxi + displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
- duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
enddo
xixl = xix(i,j,ispec)
Modified: seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -69,8 +69,9 @@
logical, dimension(nspec) :: elastic,poroelastic
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_acoustic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: veloc_elastic,velocs_poroelastic
- double precision, dimension(NDIM,npoin) :: vector_field_display
+ real(kind=CUSTOM_REAL), dimension(3,npoin) :: veloc_elastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: velocs_poroelastic
+ double precision, dimension(3,npoin) :: vector_field_display
! array with derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
@@ -80,7 +81,7 @@
integer i,j,ispec,iglob
! vector field in this element
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLX) :: vector_field_element
! loop over spectral elements
do ispec = 1,nspec
@@ -131,11 +132,12 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
! vector field in this element
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLX) :: vector_field_element
logical, dimension(nspec) :: elastic,poroelastic
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_acoustic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: veloc_elastic,velocs_poroelastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin) :: veloc_elastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: velocs_poroelastic
! array with derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
@@ -162,6 +164,7 @@
iglob = ibool(i,j,ispec)
vector_field_element(1,i,j) = veloc_elastic(1,iglob)
vector_field_element(2,i,j) = veloc_elastic(2,iglob)
+ vector_field_element(3,i,j) = veloc_elastic(3,iglob)
enddo
enddo
@@ -170,7 +173,8 @@
do i = 1,NGLLX
iglob = ibool(i,j,ispec)
vector_field_element(1,i,j) = velocs_poroelastic(1,iglob)
- vector_field_element(2,i,j) = velocs_poroelastic(2,iglob)
+ vector_field_element(2,i,j) = ZERO
+ vector_field_element(3,i,j) = velocs_poroelastic(2,iglob)
enddo
enddo
@@ -209,7 +213,8 @@
! derivatives of potential
vector_field_element(1,i,j) = (tempx1l*xixl + tempx2l*gammaxl) / rhol
- vector_field_element(2,i,j) = (tempx1l*xizl + tempx2l*gammazl) / rhol
+ vector_field_element(2,i,j) = ZERO
+ vector_field_element(3,i,j) = (tempx1l*xizl + tempx2l*gammazl) / rhol
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -279,6 +279,8 @@
double precision :: Q0,freq0
+ logical :: body_waves
+
logical, dimension(:), allocatable :: enreg_surf
integer, external :: num_4, num_9
@@ -467,6 +469,9 @@
call read_value_double_precision(IIN,IGNORE_JUNK,Q0)
call read_value_double_precision(IIN,IGNORE_JUNK,freq0)
+! determine if body or surface (membrane) waves calculation
+ call read_value_logical(IIN,IGNORE_JUNK,body_waves)
+
if ( read_external_mesh ) then
call read_mesh(mesh_file, nelmnts, elmnts, nnodes, num_start)
@@ -1444,6 +1449,9 @@
write(15,*) 'TURN_VISCATTENUATION_ON Q0 freq0'
write(15,*) TURN_VISCATTENUATION_ON,Q0,freq0
+ write(15,*) 'body_waves'
+ write(15,*) body_waves
+
write(15,*) 'nt deltat isolver'
write(15,*) nt,deltat,isolver
write(15,*) 'NSOURCE'
Modified: seismo/2D/SPECFEM2D/trunk/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotpost.F90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/plotpost.F90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -106,7 +106,7 @@
double precision dt,timeval
double precision, dimension(NSOURCE) :: x_source,z_source
- double precision displ(NDIM,npoin),coord(NDIM,npoin)
+ double precision displ(3,npoin),coord(NDIM,npoin)
double precision vpext(NGLLX,NGLLZ,nspec)
double precision coorg(NDIM,npgeo)
@@ -1464,7 +1464,7 @@
ratio_page = min(rpercentz*sizez/(zmax-zmin),rpercentx*sizex/(xmax-xmin)) / 100.d0
! compute the maximum of the norm of the vector
- dispmax = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
+ dispmax = maxval(sqrt(displ(1,:)**2 + displ(3,:)**2))
#ifdef USE_MPI
call MPI_ALLREDUCE (dispmax, dispmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
dispmax = dispmax_glob
@@ -2745,7 +2745,7 @@
do l=1,NGLLX
Uxinterp(i,j) = Uxinterp(i,j) + displ(1,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
- Uzinterp(i,j) = Uzinterp(i,j) + displ(2,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
+ Uzinterp(i,j) = Uzinterp(i,j) + displ(3,ibool(k,l,ispec))*flagrange(k,i)*flagrange(l,j)
enddo
enddo
@@ -2888,7 +2888,7 @@
z1 =(coord(2,ipoin)-zmin)*ratio_page
x2 = displ(1,ipoin)*sizemax_arrows/dispmax
- z2 = displ(2,ipoin)*sizemax_arrows/dispmax
+ z2 = displ(3,ipoin)*sizemax_arrows/dispmax
d = sqrt(x2**2 + z2**2)
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -148,7 +148,10 @@
! volume = {5336},
! pages = {350-363}}
!
-! version 6.0, Christina Morency and Yang Luo, August 2009:
+! version 6.0.1, Christina Morency, September 2009:
+! - added SH (membrane) waves calculation for elastic media
+! - fixed some bugs on acoustic kernels
+! version 6.0.0, Christina Morency and Yang Luo, August 2009:
! - support for poroelastic media
! - adjoint method for acoustic/elastic/poroelastic
!
@@ -229,6 +232,9 @@
double precision, dimension(:,:), allocatable :: coorg
double precision, dimension(:), allocatable :: coorgread
+! for body or surface (membrane) waves calculation
+ logical :: body_waves
+
! receiver information
integer :: nrec,ios
integer, dimension(:), allocatable :: ispec_selected_rec
@@ -240,7 +246,7 @@
integer :: seismo_offset, seismo_current
! vector field in an element
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLX) :: vector_field_element
! pressure in an element
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: pressure_element
@@ -250,7 +256,7 @@
integer :: i,j,k,l,it,irec,ipoin,ip,id,n,ispec,npoin,npgeo,iglob
logical :: anyabs
- double precision :: dxd,dzd,dcurld,valux,valuz,valcurl,hlagrange,rhol,cosrot,sinrot,xi,gamma,x,z
+ double precision :: dxd,dyd,dzd,dcurld,valux,valuy,valuz,valcurl,hlagrange,rhol,cosrot,sinrot,xi,gamma,x,z
! coefficients of the explicit Newmark time scheme
integer NSTEP
@@ -423,13 +429,13 @@
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_accelw_poroelastic,b_velocw_poroelastic,b_displw_poroelastic
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_accel_elastic,b_veloc_elastic,b_displ_elastic
real(kind=CUSTOM_REAL), dimension(:), allocatable :: b_potential_dot_dot_acoustic,b_potential_dot_acoustic,b_potential_acoustic
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_ac,b_displ_ac
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rho_kl, mu_kl, kappa_kl
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhol_global, mul_global, kappal_global
real(kind=CUSTOM_REAL), dimension(:), allocatable :: mu_k, kappa_k,rho_k
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhop_kl, beta_kl, alpha_kl
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rho_ac_kl, kappa_ac_kl
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhol_ac_global, kappal_ac_global
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: kappa_ac_k,rho_ac_k
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhop_ac_kl, alpha_ac_kl
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhot_kl, rhof_kl, sm_kl, eta_kl, mufr_kl, B_kl, &
C_kl, M_kl, rhob_kl, rhofb_kl, phi_kl, Bb_kl, Cb_kl, Mb_kl, mufrb_kl, &
@@ -441,7 +447,7 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: permlxx_global,permlxz_global,permlzz_global
character(len=150) :: adj_source_file
integer :: irec_local,nadj_rec_local
- double precision :: xx,zz,rholb
+ double precision :: xx,zz,rholb,tempx1l,tempx2l,b_tempx1l,b_tempx2l
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: adj_sourcearray
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_absorb_elastic_left,b_absorb_poro_s_left,b_absorb_poro_w_left
@@ -568,6 +574,7 @@
integer , dimension(:), allocatable :: left_bound,right_bound,bot_bound
double precision , dimension(:,:), allocatable :: v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot
double precision , dimension(:,:), allocatable :: t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_paco,veloc_paco,displ_paco
integer count_left,count_right,count_bot,ibegin,iend
logical :: over_critical_angle
@@ -748,6 +755,9 @@
read(IIN,"(a80)") datlin
read(IIN,*) TURN_VISCATTENUATION_ON,Q0,freq0
+ read(IIN,"(a80)") datlin
+ read(IIN,*) body_waves
+
!---- check parameters read
if (myrank == 0 .and. ipass == 1) then
write(IOUT,200) npgeo,NDIM
@@ -1013,6 +1023,20 @@
enddo !do ispec = 1,nspec
+
+ if(.not. body_waves .and. .not. any_elastic) then
+ print*, '*************** WARNING ***************'
+ print*, 'Surface (membrane) waves calculation needs an elastic medium'
+ print*, '*************** WARNING ***************'
+ stop
+ endif
+ if(body_waves .and. (TURN_ATTENUATION_ON .or. TURN_ANISOTROPY_ON)) then
+ print*, '*************** WARNING ***************'
+ print*, 'Attenuation and anisotropy are not implemented for surface (membrane) waves calculation'
+ print*, '*************** WARNING ***************'
+ stop
+ endif
+
if(TURN_ATTENUATION_ON) then
nspec_allocate = nspec
else
@@ -1170,10 +1194,10 @@
! Files to save absorbed waves needed to reconstruct backward wavefield for adjoint method
if(ipass == 1) then
if(any_elastic .and. (save_forward .or. isolver == 2)) then
- allocate(b_absorb_elastic_left(NDIM,NGLLZ,nspec_xmin,NSTEP))
- allocate(b_absorb_elastic_right(NDIM,NGLLZ,nspec_xmax,NSTEP))
- allocate(b_absorb_elastic_bottom(NDIM,NGLLX,nspec_zmin,NSTEP))
- allocate(b_absorb_elastic_top(NDIM,NGLLX,nspec_zmax,NSTEP))
+ allocate(b_absorb_elastic_left(3,NGLLZ,nspec_xmin,NSTEP))
+ allocate(b_absorb_elastic_right(3,NGLLZ,nspec_xmax,NSTEP))
+ allocate(b_absorb_elastic_bottom(3,NGLLX,nspec_zmin,NSTEP))
+ allocate(b_absorb_elastic_top(3,NGLLX,nspec_zmax,NSTEP))
endif
if(any_poroelastic .and. (save_forward .or. isolver == 2)) then
allocate(b_absorb_poro_s_left(NDIM,NGLLZ,nspec_xmin,NSTEP))
@@ -1520,7 +1544,7 @@
allocate(coord(NDIM,npoin))
! to display acoustic elements
- allocate(vector_field_display(NDIM,npoin))
+ allocate(vector_field_display(3,npoin))
if(assign_external_model) then
allocate(vpext(NGLLX,NGLLZ,nspec))
@@ -1733,8 +1757,8 @@
nadj_rec_local = nadj_rec_local + 1
endif
enddo
- if(ipass == 1) allocate(adj_sourcearray(NSTEP,NDIM,NGLLX,NGLLZ))
- if (nadj_rec_local > 0 .and. ipass == 1) allocate(adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLZ))
+ if(ipass == 1) allocate(adj_sourcearray(NSTEP,3,NGLLX,NGLLZ))
+ if (nadj_rec_local > 0 .and. ipass == 1) allocate(adj_sourcearrays(nadj_rec_local,NSTEP,3,NGLLX,NGLLZ))
irec_local = 0
do irec = 1, nrec
! compute only adjoint source arrays in the local proc
@@ -2011,9 +2035,9 @@
if(ipass == 1) then
if(any_elastic) then
- allocate(displ_elastic(NDIM,npoin))
- allocate(veloc_elastic(NDIM,npoin))
- allocate(accel_elastic(NDIM,npoin))
+ allocate(displ_elastic(3,npoin))
+ allocate(veloc_elastic(3,npoin))
+ allocate(accel_elastic(3,npoin))
allocate(rmass_inverse_elastic(npoin))
else
! allocate unused arrays with fictitious size
@@ -2024,9 +2048,9 @@
endif
! extra array if adjoint and kernels calculation
if(isolver == 2 .and. any_elastic) then
- allocate(b_displ_elastic(NDIM,npoin))
- allocate(b_veloc_elastic(NDIM,npoin))
- allocate(b_accel_elastic(NDIM,npoin))
+ allocate(b_displ_elastic(3,npoin))
+ allocate(b_veloc_elastic(3,npoin))
+ allocate(b_accel_elastic(3,npoin))
allocate(rho_kl(npoin))
allocate(rho_k(npoin))
allocate(rhol_global(npoin))
@@ -2197,11 +2221,11 @@
allocate(b_potential_acoustic(npoin))
allocate(b_potential_dot_acoustic(npoin))
allocate(b_potential_dot_dot_acoustic(npoin))
+ allocate(b_displ_ac(2,npoin))
+ allocate(accel_ac(2,npoin))
allocate(rho_ac_kl(npoin))
- allocate(rho_ac_k(npoin))
allocate(rhol_ac_global(npoin))
allocate(kappa_ac_kl(npoin))
- allocate(kappa_ac_k(npoin))
allocate(kappal_ac_global(npoin))
allocate(rhop_ac_kl(npoin))
allocate(alpha_ac_kl(npoin))
@@ -2210,11 +2234,11 @@
allocate(b_potential_acoustic(1))
allocate(b_potential_dot_acoustic(1))
allocate(b_potential_dot_dot_acoustic(1))
+ allocate(b_displ_ac(1,1))
+ allocate(accel_ac(1,1))
allocate(rho_ac_kl(1))
- allocate(rho_ac_k(1))
allocate(rhol_ac_global(1))
allocate(kappa_ac_kl(1))
- allocate(kappa_ac_k(1))
allocate(kappal_ac_global(1))
allocate(rhop_ac_kl(1))
allocate(alpha_ac_kl(1))
@@ -2994,44 +3018,92 @@
!--- left absorbing boundary
if(nspec_xmin >0) then
do ispec = 1,nspec_xmin
+
+ if(body_waves)then!P-SV waves
do id =1,2
do i=1,NGLLZ
read(35) b_absorb_elastic_left(id,i,ispec,it)
enddo
enddo
+ b_absorb_elastic_left(3,:,ispec,it) = b_absorb_elastic_left(2,:,ispec,it)
+ b_absorb_elastic_left(2,:,ispec,it) = ZERO
+ else!SH (membrane) waves
+ do i=1,NGLLZ
+ read(35) b_absorb_elastic_left(2,i,ispec,it)
+ enddo
+ b_absorb_elastic_left(1,:,ispec,it) = ZERO
+ b_absorb_elastic_left(3,:,ispec,it) = ZERO
+ endif
+
enddo
endif
!--- right absorbing boundary
if(nspec_xmax >0) then
do ispec = 1,nspec_xmax
+
+ if(body_waves)then!P-SV waves
do id =1,2
do i=1,NGLLZ
read(36) b_absorb_elastic_right(id,i,ispec,it)
enddo
enddo
+ b_absorb_elastic_right(3,:,ispec,it) = b_absorb_elastic_right(2,:,ispec,it)
+ b_absorb_elastic_right(2,:,ispec,it) = ZERO
+ else!SH (membrane) waves
+ do i=1,NGLLZ
+ read(36) b_absorb_elastic_right(2,i,ispec,it)
+ enddo
+ b_absorb_elastic_right(1,:,ispec,it) = ZERO
+ b_absorb_elastic_right(3,:,ispec,it) = ZERO
+ endif
+
enddo
endif
!--- bottom absorbing boundary
if(nspec_zmin >0) then
do ispec = 1,nspec_zmin
+
+ if(body_waves)then!P-SV waves
do id =1,2
do i=1,NGLLX
read(37) b_absorb_elastic_bottom(id,i,ispec,it)
enddo
enddo
+ b_absorb_elastic_bottom(3,:,ispec,it) = b_absorb_elastic_bottom(2,:,ispec,it)
+ b_absorb_elastic_bottom(2,:,ispec,it) = ZERO
+ else!SH (membrane) waves
+ do i=1,NGLLZ
+ read(37) b_absorb_elastic_bottom(2,i,ispec,it)
+ enddo
+ b_absorb_elastic_bottom(1,:,ispec,it) = ZERO
+ b_absorb_elastic_bottom(3,:,ispec,it) = ZERO
+ endif
+
enddo
endif
!--- top absorbing boundary
if(nspec_zmax >0) then
do ispec = 1,nspec_zmax
+
+ if(body_waves)then!P-SV waves
do id =1,2
do i=1,NGLLX
read(38) b_absorb_elastic_top(id,i,ispec,it)
enddo
enddo
+ b_absorb_elastic_top(3,:,ispec,it) = b_absorb_elastic_top(2,:,ispec,it)
+ b_absorb_elastic_top(2,:,ispec,it) = ZERO
+ else!SH (membrane) waves
+ do i=1,NGLLZ
+ read(38) b_absorb_elastic_top(2,i,ispec,it)
+ enddo
+ b_absorb_elastic_top(1,:,ispec,it) = ZERO
+ b_absorb_elastic_top(3,:,ispec,it) = ZERO
+ endif
+
enddo
endif
@@ -3156,11 +3228,31 @@
if(any_elastic) then
write(outputname,'(a,i6.6,a)') 'lastframe_elastic',myrank,'.bin'
open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
+ if(body_waves)then !P-SV waves
do j=1,npoin
read(55) (b_displ_elastic(i,j), i=1,NDIM), &
(b_veloc_elastic(i,j), i=1,NDIM), &
(b_accel_elastic(i,j), i=1,NDIM)
enddo
+ b_displ_elastic(3,:) = b_displ_elastic(2,:)
+ b_displ_elastic(2,:) = ZERO
+ b_veloc_elastic(3,:) = b_veloc_elastic(2,:)
+ b_veloc_elastic(2,:) = ZERO
+ b_accel_elastic(3,:) = b_accel_elastic(2,:)
+ b_accel_elastic(2,:) = ZERO
+ else !SH (membrane) waves
+ do j=1,npoin
+ read(55) b_displ_elastic(2,j), &
+ b_veloc_elastic(2,j), &
+ b_accel_elastic(2,j)
+ enddo
+ b_displ_elastic(1,:) = ZERO
+ b_displ_elastic(3,:) = ZERO
+ b_veloc_elastic(1,:) = ZERO
+ b_veloc_elastic(3,:) = ZERO
+ b_accel_elastic(1,:) = ZERO
+ b_accel_elastic(3,:) = ZERO
+ endif
close(55)
write(outputname,'(a,i6.6,a)') 'snapshot_rho_kappa_mu_',myrank
@@ -3468,7 +3560,7 @@
displ_elastic(1,i) = A_plane(1) * ricker_Bielak_displ(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ B_plane(1) * ricker_Bielak_displ(t - sin(angleforce(1))*x/c_inc - cos(angleforce(1))*z/c_inc,f0(1)) &
+ C_plane(1) * ricker_Bielak_displ(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0(1))
- displ_elastic(2,i) = A_plane(2) * ricker_Bielak_displ(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ displ_elastic(3,i) = A_plane(2) * ricker_Bielak_displ(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ B_plane(2) * ricker_Bielak_displ(t - sin(angleforce(1))*x/c_inc - cos(angleforce(1))*z/c_inc,f0(1)) &
+ C_plane(2) * ricker_Bielak_displ(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0(1))
@@ -3476,7 +3568,7 @@
veloc_elastic(1,i) = A_plane(1) * ricker_Bielak_veloc(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ B_plane(1) * ricker_Bielak_veloc(t - sin(angleforce(1))*x/c_inc - cos(angleforce(1))*z/c_inc,f0(1)) &
+ C_plane(1) * ricker_Bielak_veloc(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0(1))
- veloc_elastic(2,i) = A_plane(2) * ricker_Bielak_veloc(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ veloc_elastic(3,i) = A_plane(2) * ricker_Bielak_veloc(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ B_plane(2) * ricker_Bielak_veloc(t - sin(angleforce(1))*x/c_inc - cos(angleforce(1))*z/c_inc,f0(1)) &
+ C_plane(2) * ricker_Bielak_veloc(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0(1))
@@ -3484,7 +3576,7 @@
accel_elastic(1,i) = A_plane(1) * ricker_Bielak_accel(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ B_plane(1) * ricker_Bielak_accel(t - sin(angleforce(1))*x/c_inc - cos(angleforce(1))*z/c_inc,f0(1)) &
+ C_plane(1) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0(1))
- accel_elastic(2,i) = A_plane(2) * ricker_Bielak_accel(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ accel_elastic(3,i) = A_plane(2) * ricker_Bielak_accel(t - sin(angleforce(1))*x/c_inc + cos(angleforce(1))*z/c_inc,f0(1)) &
+ B_plane(2) * ricker_Bielak_accel(t - sin(angleforce(1))*x/c_inc - cos(angleforce(1))*z/c_inc,f0(1)) &
+ C_plane(2) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0(1))
@@ -3556,17 +3648,32 @@
allocate(t0x_bot(count_bot,NSTEP))
allocate(t0z_bot(count_bot,NSTEP))
+ allocate(displ_paco(NDIM,npoin))
+ allocate(veloc_paco(NDIM,npoin))
+ allocate(accel_paco(NDIM,npoin))
+
! call Paco's routine to compute in frequency and convert to time by Fourier transform
call paco_beyond_critical(coord,npoin,deltat,NSTEP,angleforce(1),&
f0(1),cploc,csloc,TURN_ATTENUATION_ON,Qp_attenuation,source_type(1),v0x_left,v0z_left,&
v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,&
t0x_bot,t0z_bot,left_bound(1:count_left),right_bound(1:count_right),bot_bound(1:count_bot)&
- ,count_left,count_right,count_bot,displ_elastic,veloc_elastic,accel_elastic)
+ ,count_left,count_right,count_bot,displ_paco,veloc_paco,accel_paco)
+ displ_elastic(1,:) = displ_paco(1,:)
+ displ_elastic(3,:) = displ_paco(2,:)
+ veloc_elastic(1,:) = veloc_paco(1,:)
+ veloc_elastic(3,:) = veloc_paco(2,:)
+ accel_elastic(1,:) = accel_paco(1,:)
+ accel_elastic(3,:) = accel_paco(2,:)
+
deallocate(left_bound)
deallocate(right_bound)
deallocate(bot_bound)
+ deallocate(displ_paco)
+ deallocate(veloc_paco)
+ deallocate(accel_paco)
+
if (myrank == 0) then
write(IOUT,*) '***********'
write(IOUT,*) 'done calculating the initial wave field'
@@ -3575,7 +3682,7 @@
endif ! beyond critical angle
- write(IOUT,*) 'Max norm of initial elastic displacement = ',maxval(sqrt(displ_elastic(1,:)**2 + displ_elastic(2,:)**2))
+ write(IOUT,*) 'Max norm of initial elastic displacement = ',maxval(sqrt(displ_elastic(1,:)**2 + displ_elastic(3,:)**2))
endif ! initialfield
deltatsquare = deltat * deltat
@@ -3608,10 +3715,10 @@
! Ricker (second derivative of a Gaussian) source time function
if(time_function_type(i_source) == 1) then
- source_time_function(i_source,it) = - factor(i_source) * (ONE-TWO*aval(i_source)*(time-t0(i_source))**2) * &
- exp(-aval(i_source)*(time-t0(i_source))**2)
-! source_time_function(i_source,it) = - factor(i_source) * TWO*aval(i_source)*sqrt(aval(i_source))*&
-! (time-t0(i_source))/pi * exp(-aval(i_source)*(time-t0(i_source))**2)
+! source_time_function(i_source,it) = - factor(i_source) * (ONE-TWO*aval(i_source)*(time-t0(i_source))**2) * &
+! exp(-aval(i_source)*(time-t0(i_source))**2)
+ source_time_function(i_source,it) = - factor(i_source) * TWO*aval(i_source)*sqrt(aval(i_source))*&
+ (time-t0(i_source))/pi * exp(-aval(i_source)*(time-t0(i_source))**2)
! first derivative of a Gaussian source time function
else if(time_function_type(i_source) == 2) then
@@ -4687,7 +4794,7 @@
jbegin_left,jend_left,jbegin_right,jend_right,isolver,save_forward,b_absorb_acoustic_left,&
b_absorb_acoustic_right,b_absorb_acoustic_bottom,&
b_absorb_acoustic_top,nspec_xmin,nspec_xmax,&
- nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax,kappa_ac_k)
+ nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax)
if(anyabs .and. save_forward .and. isolver == 1) then
@@ -4757,11 +4864,11 @@
iglob = ibool(i,j,ispec_elastic)
displ_x = displ_elastic(1,iglob)
- displ_z = displ_elastic(2,iglob)
+ displ_z = displ_elastic(3,iglob)
if(isolver == 2) then
b_displ_x = b_displ_elastic(1,iglob)
- b_displ_z = b_displ_elastic(2,iglob)
+ b_displ_z = b_displ_elastic(3,iglob)
endif
! get point values for the acoustic side
@@ -5009,19 +5116,13 @@
endif
endif
- if(any_acoustic .and. isolver == 2) then ! kernels calculation
- do iglob = 1,npoin
- rho_ac_k(iglob) = potential_dot_dot_acoustic(iglob)*b_potential_acoustic(iglob)
- enddo
- endif
-
! *********************************************************
! ************* main solver for the elastic elements
! *********************************************************
if(any_elastic) then
- call compute_forces_elastic(npoin,nspec,myrank,nelemabs,numat, &
+ call compute_forces_elastic(body_waves,npoin,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
@@ -5044,44 +5145,81 @@
!--- left absorbing boundary
if(nspec_xmin >0) then
do ispec = 1,nspec_xmin
- do id =1,2
+
+ if(body_waves)then!P-SV waves
do i=1,NGLLZ
- write(35) b_absorb_elastic_left(id,i,ispec,it)
+ write(35) b_absorb_elastic_left(1,i,ispec,it)
enddo
- enddo
+ do i=1,NGLLZ
+ write(35) b_absorb_elastic_left(3,i,ispec,it)
+ enddo
+ else!SH (membrane) waves
+ do i=1,NGLLZ
+ write(35) b_absorb_elastic_left(2,i,ispec,it)
+ enddo
+ endif
+
enddo
endif
!--- right absorbing boundary
if(nspec_xmax >0) then
do ispec = 1,nspec_xmax
- do id =1,2
+
+
+ if(body_waves)then!P-SV waves
do i=1,NGLLZ
- write(36) b_absorb_elastic_right(id,i,ispec,it)
+ write(36) b_absorb_elastic_right(1,i,ispec,it)
enddo
- enddo
+ do i=1,NGLLZ
+ write(36) b_absorb_elastic_right(3,i,ispec,it)
+ enddo
+ else!SH (membrane) waves
+ do i=1,NGLLZ
+ write(36) b_absorb_elastic_right(2,i,ispec,it)
+ enddo
+ endif
+
enddo
endif
!--- bottom absorbing boundary
if(nspec_zmin >0) then
do ispec = 1,nspec_zmin
- do id =1,2
+
+ if(body_waves)then!P-SV waves
do i=1,NGLLX
- write(37) b_absorb_elastic_bottom(id,i,ispec,it)
+ write(37) b_absorb_elastic_bottom(1,i,ispec,it)
enddo
- enddo
+ do i=1,NGLLX
+ write(37) b_absorb_elastic_bottom(3,i,ispec,it)
+ enddo
+ else!SH (membrane) waves
+ do i=1,NGLLX
+ write(37) b_absorb_elastic_bottom(2,i,ispec,it)
+ enddo
+ endif
+
enddo
endif
!--- top absorbing boundary
if(nspec_zmax >0) then
do ispec = 1,nspec_zmax
- do id =1,2
+
+ if(body_waves)then!P-SV waves
do i=1,NGLLX
- write(38) b_absorb_elastic_top(id,i,ispec,it)
+ write(38) b_absorb_elastic_top(1,i,ispec,it)
enddo
- enddo
+ do i=1,NGLLX
+ write(38) b_absorb_elastic_top(3,i,ispec,it)
+ enddo
+ else!SH (membrane) waves
+ do i=1,NGLLX
+ write(38) b_absorb_elastic_top(2,i,ispec,it)
+ enddo
+ endif
+
enddo
endif
@@ -5160,11 +5298,11 @@
endif
accel_elastic(1,iglob) = accel_elastic(1,iglob) + weight*nx*pressure
- accel_elastic(2,iglob) = accel_elastic(2,iglob) + weight*nz*pressure
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + weight*nz*pressure
if(isolver == 2) then
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + weight*nx*b_pressure
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) + weight*nz*b_pressure
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) + weight*nz*b_pressure
endif !if(isolver == 2) then
enddo
@@ -5348,15 +5486,15 @@
! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
dux_dxi = dux_dxi + displ_elastic(1,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
- duz_dxi = duz_dxi + displ_elastic(2,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
+ duz_dxi = duz_dxi + displ_elastic(3,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
dux_dgamma = dux_dgamma + displ_elastic(1,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
- duz_dgamma = duz_dgamma + displ_elastic(2,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
+ duz_dgamma = duz_dgamma + displ_elastic(3,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
if(isolver == 2) then
b_dux_dxi = b_dux_dxi + b_displ_elastic(1,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
- b_duz_dxi = b_duz_dxi + b_displ_elastic(2,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
+ b_duz_dxi = b_duz_dxi + b_displ_elastic(3,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
b_dux_dgamma = b_dux_dgamma + b_displ_elastic(1,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
- b_duz_dgamma = b_duz_dgamma + b_displ_elastic(2,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
+ b_duz_dgamma = b_duz_dgamma + b_displ_elastic(3,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
endif
enddo
@@ -5437,14 +5575,14 @@
accel_elastic(1,iglob) = accel_elastic(1,iglob) - weight* &
(sigma_xx*nx + sigma_xz*nz +sigmap*nx)/3.d0
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - weight* &
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - weight* &
(sigma_xz*nx + sigma_zz*nz +sigmap*nz)/3.d0
if(isolver == 2) then
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - weight* &
(b_sigma_xx*nx + b_sigma_xz*nz +b_sigmap*nx)/3.d0
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - weight* &
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - weight* &
(b_sigma_xz*nx + b_sigma_zz*nz +b_sigmap*nz)/3.d0
endif !if(isolver == 2) then
@@ -5493,15 +5631,29 @@
! collocated force
if(source_type(i_source) == 1) then
if(isolver == 1) then ! forward wavefield
+
+ if(body_waves) then ! P-SV calculation
accel_elastic(1,iglob_source(i_source)) = accel_elastic(1,iglob_source(i_source)) &
- sin(angleforce(i_source))*source_time_function(i_source,it)
+ accel_elastic(3,iglob_source(i_source)) = accel_elastic(3,iglob_source(i_source)) &
+ + cos(angleforce(i_source))*source_time_function(i_source,it)
+ else ! SH (membrane) calculation
accel_elastic(2,iglob_source(i_source)) = accel_elastic(2,iglob_source(i_source)) &
- + cos(angleforce(i_source))*source_time_function(i_source,it)
+ + source_time_function(i_source,it)
+ endif
+
else ! backward wavefield
+
+ if(body_waves) then ! P-SV calculation
b_accel_elastic(1,iglob_source(i_source)) = b_accel_elastic(1,iglob_source(i_source)) &
- sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ b_accel_elastic(3,iglob_source(i_source)) = b_accel_elastic(3,iglob_source(i_source)) &
+ + cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ else ! SH (membrane) calculation
b_accel_elastic(2,iglob_source(i_source)) = b_accel_elastic(2,iglob_source(i_source)) &
- + cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ + source_time_function(i_source,NSTEP-it+1)
+ endif
+
endif !endif isolver == 1
endif
@@ -5512,24 +5664,20 @@
accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
+ accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic
veloc_elastic = veloc_elastic + deltatover2*accel_elastic
if(isolver == 2) then
b_accel_elastic(1,:) = b_accel_elastic(1,:) * rmass_inverse_elastic(:)
b_accel_elastic(2,:) = b_accel_elastic(2,:) * rmass_inverse_elastic(:)
+ b_accel_elastic(3,:) = b_accel_elastic(3,:) * rmass_inverse_elastic(:)
b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic
endif
endif
- if(any_elastic .and. isolver == 2) then ! kernels calculation
- do iglob = 1,npoin
- rho_k(iglob) = accel_elastic(1,iglob)*b_displ_elastic(1,iglob) +&
- accel_elastic(2,iglob)*b_displ_elastic(2,iglob)
- enddo
- endif
! ******************************************************************************************************************
! ************* main solver for the poroelastic elements: first the solid (u_s) than the fluid (w)
@@ -5801,15 +5949,15 @@
! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec_elastic))*hprime_xx(i,k)
- duz_dxi = duz_dxi + displ_elastic(2,ibool(k,j,ispec_elastic))*hprime_xx(i,k)
+ duz_dxi = duz_dxi + displ_elastic(3,ibool(k,j,ispec_elastic))*hprime_xx(i,k)
dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec_elastic))*hprime_zz(j,k)
- duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec_elastic))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displ_elastic(3,ibool(i,k,ispec_elastic))*hprime_zz(j,k)
if(isolver == 2) then
b_dux_dxi = b_dux_dxi + b_displ_elastic(1,ibool(k,j,ispec_elastic))*hprime_xx(i,k)
- b_duz_dxi = b_duz_dxi + b_displ_elastic(2,ibool(k,j,ispec_elastic))*hprime_xx(i,k)
+ b_duz_dxi = b_duz_dxi + b_displ_elastic(3,ibool(k,j,ispec_elastic))*hprime_xx(i,k)
b_dux_dgamma = b_dux_dgamma + b_displ_elastic(1,ibool(i,k,ispec_elastic))*hprime_zz(j,k)
- b_duz_dgamma = b_duz_dgamma + b_displ_elastic(2,ibool(i,k,ispec_elastic))*hprime_zz(j,k)
+ b_duz_dgamma = b_duz_dgamma + b_displ_elastic(3,ibool(i,k,ispec_elastic))*hprime_zz(j,k)
endif
enddo
@@ -6142,20 +6290,6 @@
endif
- if(any_poroelastic .and. isolver ==2) then
- do iglob =1,npoin
- rhot_k(iglob) = accels_poroelastic(1,iglob) * b_displs_poroelastic(1,iglob) + &
- accels_poroelastic(2,iglob) * b_displs_poroelastic(2,iglob)
- rhof_k(iglob) = accelw_poroelastic(1,iglob) * b_displs_poroelastic(1,iglob) + &
- accelw_poroelastic(2,iglob) * b_displs_poroelastic(2,iglob) + &
- accels_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
- accels_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
- sm_k(iglob) = accelw_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
- accelw_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
- eta_k(iglob) = velocw_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
- velocw_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
- enddo
- endif
!assembling the displacements on the elastic-poro boundaries
if(coupled_elastic_poro) then
@@ -6179,9 +6313,9 @@
if(icount(iglob) ==1)then
veloc_elastic(1,iglob) = veloc_elastic(1,iglob) - deltatover2*accel_elastic(1,iglob)
- veloc_elastic(2,iglob) = veloc_elastic(2,iglob) - deltatover2*accel_elastic(2,iglob)
+ veloc_elastic(3,iglob) = veloc_elastic(3,iglob) - deltatover2*accel_elastic(3,iglob)
accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
- accel_elastic(2,iglob) = accel_elastic(2,iglob) / rmass_inverse_elastic(iglob)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
! recovering original velocities and accelerations on boundaries (poro side)
velocs_poroelastic(1,iglob) = velocs_poroelastic(1,iglob) - deltatover2*accels_poroelastic(1,iglob)
velocs_poroelastic(2,iglob) = velocs_poroelastic(2,iglob) - deltatover2*accels_poroelastic(2,iglob)
@@ -6190,15 +6324,15 @@
! assembling accelerations
accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
- accel_elastic(2,iglob) = ( accel_elastic(2,iglob) + accels_poroelastic(2,iglob) ) / &
+ accel_elastic(3,iglob) = ( accel_elastic(3,iglob) + accels_poroelastic(2,iglob) ) / &
( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
- accels_poroelastic(2,iglob) = accel_elastic(2,iglob)
+ accels_poroelastic(2,iglob) = accel_elastic(3,iglob)
! updating velocities
velocs_poroelastic(1,iglob) = velocs_poroelastic(1,iglob) + deltatover2*accels_poroelastic(1,iglob)
velocs_poroelastic(2,iglob) = velocs_poroelastic(2,iglob) + deltatover2*accels_poroelastic(2,iglob)
veloc_elastic(1,iglob) = veloc_elastic(1,iglob) + deltatover2*accel_elastic(1,iglob)
- veloc_elastic(2,iglob) = veloc_elastic(2,iglob) + deltatover2*accel_elastic(2,iglob)
+ veloc_elastic(3,iglob) = veloc_elastic(3,iglob) + deltatover2*accel_elastic(3,iglob)
! zeros w
accelw_poroelastic(1,iglob) = ZERO
accelw_poroelastic(2,iglob) = ZERO
@@ -6211,7 +6345,30 @@
enddo
endif
+! kernels calculation
+ if(any_elastic .and. isolver == 2) then ! kernels calculation
+ do iglob = 1,npoin
+ rho_k(iglob) = accel_elastic(1,iglob)*b_displ_elastic(1,iglob) +&
+ accel_elastic(2,iglob)*b_displ_elastic(2,iglob) +&
+ accel_elastic(3,iglob)*b_displ_elastic(3,iglob)
+ enddo
+ endif
+ if(any_poroelastic .and. isolver ==2) then
+ do iglob =1,npoin
+ rhot_k(iglob) = accels_poroelastic(1,iglob) * b_displs_poroelastic(1,iglob) + &
+ accels_poroelastic(2,iglob) * b_displs_poroelastic(2,iglob)
+ rhof_k(iglob) = accelw_poroelastic(1,iglob) * b_displs_poroelastic(1,iglob) + &
+ accelw_poroelastic(2,iglob) * b_displs_poroelastic(2,iglob) + &
+ accels_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
+ accels_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
+ sm_k(iglob) = accelw_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
+ accelw_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
+ eta_k(iglob) = velocw_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
+ velocw_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
+ enddo
+ endif
+
!---- compute kinetic and potential energy
if(OUTPUT_ENERGY) &
call compute_energy(displ_elastic,veloc_elastic,displs_poroelastic,velocs_poroelastic, &
@@ -6238,7 +6395,7 @@
if(any_elastic_glob) then
if(any_elastic) then
- displnorm_all = maxval(sqrt(displ_elastic(1,:)**2 + displ_elastic(2,:)**2))
+ displnorm_all = maxval(sqrt(displ_elastic(1,:)**2 + displ_elastic(2,:)**2 + displ_elastic(3,:)**2))
else
displnorm_all = 0.d0
endif
@@ -6348,6 +6505,7 @@
! perform the general interpolation using Lagrange polynomials
valux = ZERO
+ valuy = ZERO
valuz = ZERO
valcurl = ZERO
@@ -6382,7 +6540,8 @@
dzd = displs_poroelastic(2,iglob)
elseif(elastic(ispec)) then
dxd = displ_elastic(1,iglob)
- dzd = displ_elastic(2,iglob)
+ dyd = displ_elastic(2,iglob)
+ dzd = displ_elastic(3,iglob)
endif
else if(seismotype == 2) then
@@ -6392,7 +6551,8 @@
dzd = velocs_poroelastic(2,iglob)
elseif(elastic(ispec)) then
dxd = veloc_elastic(1,iglob)
- dzd = veloc_elastic(2,iglob)
+ dyd = veloc_elastic(2,iglob)
+ dzd = veloc_elastic(3,iglob)
endif
else if(seismotype == 3) then
@@ -6402,7 +6562,8 @@
dzd = accels_poroelastic(2,iglob)
elseif(elastic(ispec)) then
dxd = accel_elastic(1,iglob)
- dzd = accel_elastic(2,iglob)
+ dyd = accel_elastic(2,iglob)
+ dzd = accel_elastic(3,iglob)
endif
else if(seismotype == 5) then
@@ -6420,6 +6581,7 @@
! compute interpolated field
valux = valux + dxd*hlagrange
+ valuy = valuy + dyd*hlagrange
valuz = valuz + dzd*hlagrange
valcurl = valcurl + dcurld*hlagrange
@@ -6428,8 +6590,13 @@
! rotate seismogram components if needed, except if recording pressure, which is a scalar
if(seismotype /= 4 .and. seismotype /= 6) then
+ if(body_waves) then
sisux(seismo_current,irecloc) = cosrot_irec(irecloc)*valux + sinrot_irec(irecloc)*valuz
sisuz(seismo_current,irecloc) = - sinrot_irec(irecloc)*valux + cosrot_irec(irecloc)*valuz
+ else
+ sisux(seismo_current,irecloc) = valuy
+ sisuz(seismo_current,irecloc) = ZERO
+ endif
else
sisux(seismo_current,irecloc) = valux
sisuz(seismo_current,irecloc) = ZERO
@@ -6449,19 +6616,52 @@
do ispec = 1, nspec
if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- do k = 1, NGLLZ
+ do j = 1, NGLLZ
do i = 1, NGLLX
- iglob = ibool(i,k,ispec)
+ iglob = ibool(i,j,ispec)
kappal_ac_global(iglob) = poroelastcoef(3,1,kmato(ispec))
rhol_ac_global(iglob) = density(1,kmato(ispec))
- enddo
- enddo
+
+! calcul the displacement by computing the gradient of potential / rho
+! and calcul the acceleration by computing the gradient of potential_dot_dot / rho
+ tempx1l = ZERO
+ tempx2l = ZERO
+ b_tempx1l = ZERO
+ b_tempx2l = ZERO
+ do k = 1,NGLLX
+! derivative along x
+ tempx1l = tempx1l + potential_dot_dot_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
+ b_tempx1l = b_tempx1l + b_potential_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
+! derivative along z
+ tempx2l = tempx2l + potential_dot_dot_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
+ b_tempx2l = b_tempx2l + b_potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
+ enddo
+
+ xixl = xix(i,j,ispec)
+ xizl = xiz(i,j,ispec)
+ gammaxl = gammax(i,j,ispec)
+ gammazl = gammaz(i,j,ispec)
+
+ if(assign_external_model) rhol_ac_global(iglob) = rhoext(i,j,ispec)
+
+! derivatives of potential
+ accel_ac(1,iglob) = (tempx1l*xixl + tempx2l*gammaxl) / rhol_ac_global(iglob)
+ accel_ac(2,iglob) = (tempx1l*xizl + tempx2l*gammazl) / rhol_ac_global(iglob)
+ b_displ_ac(1,iglob) = (b_tempx1l*xixl + b_tempx2l*gammaxl) / rhol_ac_global(iglob)
+ b_displ_ac(2,iglob) = (b_tempx1l*xizl + b_tempx2l*gammazl) / rhol_ac_global(iglob)
+
+ enddo !i = 1, NGLLX
+ enddo !j = 1, NGLLZ
endif
enddo
do iglob =1,npoin
- rho_ac_kl(iglob) = rho_ac_kl(iglob) - rhol_ac_global(iglob) * rho_ac_k(iglob) * deltat
- kappa_ac_kl(iglob) = kappa_ac_kl(iglob) - kappal_ac_global(iglob) * kappa_ac_k(iglob) * deltat
+ rho_ac_kl(iglob) = rho_ac_kl(iglob) - rhol_ac_global(iglob) * &
+ dot_product(accel_ac(:,iglob),b_displ_ac(:,iglob)) * deltat
+ kappa_ac_kl(iglob) = kappa_ac_kl(iglob) - kappal_ac_global(iglob) * &
+ potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob) * &
+ b_potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob)&
+ * deltat
!
rhop_ac_kl(iglob) = rho_ac_kl(iglob) + kappa_ac_kl(iglob)
alpha_ac_kl(iglob) = TWO * kappa_ac_kl(iglob)
@@ -6715,7 +6915,7 @@
if (myrank == 0) write(IOUT,*) 'Writing PostScript file'
- if(imagetype == 1) then
+ if(imagetype == 1 .and. body_waves) then
if (myrank == 0) write(IOUT,*) 'drawing displacement vector as small arrows...'
@@ -6751,7 +6951,7 @@
d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field,d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
coorg_send_ps_vector_field,coorg_recv_ps_vector_field)
- else if(imagetype == 2) then
+ else if(imagetype == 2 .and. body_waves) then
if (myrank == 0) write(IOUT,*) 'drawing velocity vector as small arrows...'
@@ -6787,7 +6987,7 @@
d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field,d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
coorg_send_ps_vector_field,coorg_recv_ps_vector_field)
- else if(imagetype == 3) then
+ else if(imagetype == 3 .and. body_waves) then
if (myrank == 0) write(IOUT,*) 'drawing acceleration vector as small arrows...'
@@ -6823,15 +7023,15 @@
d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field,d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
coorg_send_ps_vector_field,coorg_recv_ps_vector_field)
- else if(imagetype == 4) then
+ else if(imagetype == 4 .or. .not. body_waves) then
- if (myrank == 0) write(IOUT,*) 'cannot draw scalar pressure field as a vector plot, skipping...'
+ if (myrank == 0) write(IOUT,*) 'cannot draw scalar pressure field or y-component field as a vector plot, skipping...'
else
call exit_MPI('wrong type for snapshots')
endif
- if (myrank == 0 .and. imagetype /= 4) write(IOUT,*) 'PostScript file written'
+ if (myrank == 0 .and. imagetype /= 4 .and. body_waves) write(IOUT,*) 'PostScript file written'
endif
@@ -6844,7 +7044,7 @@
if(imagetype == 1) then
- if (myrank == 0) write(IOUT,*) 'drawing image of vertical component of displacement vector...'
+ if (myrank == 0) write(IOUT,*) 'drawing image of z (if P-SV) or y (if SH) component of displacement vector...'
call compute_vector_whole_medium(potential_acoustic,displ_elastic,displs_poroelastic,&
elastic,poroelastic,vector_field_display, &
@@ -6852,7 +7052,7 @@
else if(imagetype == 2) then
- if (myrank == 0) write(IOUT,*) 'drawing image of vertical component of velocity vector...'
+ if (myrank == 0) write(IOUT,*) 'drawing image of z (if P-SV) or y (if SH) component of velocity vector...'
call compute_vector_whole_medium(potential_dot_acoustic,veloc_elastic,velocs_poroelastic,&
elastic,poroelastic,vector_field_display, &
@@ -6860,13 +7060,13 @@
else if(imagetype == 3) then
- if (myrank == 0) write(IOUT,*) 'drawing image of vertical component of acceleration vector...'
+ if (myrank == 0) write(IOUT,*) 'drawing image of z (if P-SV) or y (if SH) component of acceleration vector...'
call compute_vector_whole_medium(potential_dot_dot_acoustic,accel_elastic,accels_poroelastic,&
elastic,poroelastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
- else if(imagetype == 4) then
+ else if(imagetype == 4 .and. body_waves) then
if (myrank == 0) write(IOUT,*) 'drawing image of pressure field...'
@@ -6876,6 +7076,8 @@
numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext,e1,e11, &
TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
+ else if(imagetype == 4 .and. .not. body_waves) then
+ call exit_MPI('cannot draw pressure field for SH (membrane) waves')
else
call exit_MPI('wrong type for snapshots')
endif
@@ -6885,7 +7087,11 @@
do k = 1, nb_pixel_loc
j = ceiling(real(num_pixel_loc(k)) / real(NX_IMAGE_color))
i = num_pixel_loc(k) - (j-1)*NX_IMAGE_color
+ if(body_waves) then !P-SH waves, plot vertical component or pressure
+ image_color_data(i,j) = vector_field_display(3,iglob_image_color(i,j))
+ else !SH (membrane) waves, plot y-component
image_color_data(i,j) = vector_field_display(2,iglob_image_color(i,j))
+ endif
enddo
! assembling array image_color_data on process zero for color output
@@ -6908,7 +7114,11 @@
do k = 1, nb_pixel_loc
j = ceiling(real(num_pixel_loc(k)) / real(NX_IMAGE_color))
i = num_pixel_loc(k) - (j-1)*NX_IMAGE_color
+ if(body_waves) then !P-SH waves, plot vertical component or pressure
+ data_pixel_send(k) = vector_field_display(3,iglob_image_color(i,j))
+ else !SH (membrane) waves, plot y-component
data_pixel_send(k) = vector_field_display(2,iglob_image_color(i,j))
+ endif
enddo
call MPI_SEND(data_pixel_send(1),nb_pixel_loc,MPI_DOUBLE_PRECISION, 0, 43, MPI_COMM_WORLD, ier)
@@ -6929,7 +7139,7 @@
! suppress seismograms if we generate traces of the run for analysis with "ParaVer", because time consuming
if(.not. GENERATE_PARAVER_TRACES) call write_seismograms(sisux,sisuz,siscurl,station_name,network_name,NSTEP, &
nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0(1), &
- NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current)
+ NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,body_waves)
seismo_offset = seismo_offset + seismo_current
seismo_current = 0
@@ -7003,11 +7213,19 @@
endif
write(outputname,'(a,i6.6,a)') 'lastframe_elastic',myrank,'.bin'
open(unit=55,file='OUTPUT_FILES/'//outputname,status='unknown',form='unformatted')
+ if(body_waves)then !P-SV waves
do j=1,npoin
- write(55) (displ_elastic(i,j), i=1,NDIM), &
- (veloc_elastic(i,j), i=1,NDIM), &
- (accel_elastic(i,j), i=1,NDIM)
+ write(55) displ_elastic(1,j), displ_elastic(3,j), &
+ veloc_elastic(1,j), veloc_elastic(3,j), &
+ accel_elastic(1,j), accel_elastic(3,j)
enddo
+ else !SH (membrane) waves
+ do j=1,npoin
+ write(55) displ_elastic(2,j), &
+ veloc_elastic(2,j), &
+ accel_elastic(2,j)
+ enddo
+ endif
close(55)
endif
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2009-09-11 14:15:49 UTC (rev 15665)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2009-09-14 21:56:41 UTC (rev 15666)
@@ -46,7 +46,7 @@
subroutine write_seismograms(sisux,sisuz,siscurl,station_name,network_name, &
NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
- NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current &
+ NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,body_waves &
)
implicit none
@@ -60,6 +60,8 @@
integer :: NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current
double precision :: t0,deltat
+ logical :: body_waves
+
integer, intent(in) :: nrecloc,myrank
integer, dimension(nrec),intent(in) :: which_proc_receiver
@@ -110,8 +112,8 @@
endif
-! only one seismogram if pressurs
- if(seismotype == 4 .or. seismotype == 6) then
+! only one seismogram if pressures or SH (membrane) waves
+ if(seismotype == 4 .or. seismotype == 6 .or. .not. body_waves) then
number_of_components = 1
else if(seismotype == 5) then
number_of_components = NDIM+1
@@ -156,18 +158,22 @@
! write the new files
if(seismotype == 4 .or. seismotype == 6) then
open(unit=12,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown',access='direct',recl=4)
+ elseif(.not.body_waves) then
+ open(unit=12,file='OUTPUT_FILES/Uy_file_single.bin',status='unknown',access='direct',recl=4)
else
open(unit=12,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown',access='direct',recl=4)
endif
if(seismotype == 4 .or. seismotype == 6) then
open(unit=13,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown',access='direct',recl=8)
+ elseif(.not.body_waves) then
+ open(unit=13,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown',access='direct',recl=8)
else
open(unit=13,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown',access='direct',recl=8)
endif
! no Z component seismogram if pressure
- if(seismotype /= 4 .and. seismotype /= 6) then
+ if(seismotype /= 4 .and. seismotype /= 6 .and. body_waves) then
open(unit=14,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown',access='direct',recl=4)
open(unit=15,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown',access='direct',recl=8)
@@ -232,6 +238,8 @@
! in case of pressure, use different abbreviation
if(seismotype == 4 .or. seismotype == 6) chn = 'PRE'
+ ! in case of SH (membrane) waves, use different abbreviation
+ if(.not.body_waves) chn = 'BHY'
! create the name of the seismogram file for each slice
! file name includes the name of the station, the network and the component
@@ -277,7 +285,7 @@
do isample = 1, seismo_current
write(12,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1))
write(13,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,1)
- if ( seismotype /= 4 .and. seismotype /= 6) then
+ if ( seismotype /= 4 .and. seismotype /= 6 .and. body_waves) then
write(14,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
write(15,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,2)
end if
@@ -308,7 +316,7 @@
close(12)
close(13)
- if ( seismotype /= 4 .and. seismotype /= 6) then
+ if ( seismotype /= 4 .and. seismotype /= 6 .and. body_waves) then
close(14)
close(15)
end if
More information about the CIG-COMMITS
mailing list