[cig-commits] r20342 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES/Tape2007 setup src/specfem2D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Sun Jun 10 14:11:47 PDT 2012
Author: danielpeter
Date: 2012-06-10 14:11:46 -0700 (Sun, 10 Jun 2012)
New Revision: 20342
Added:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
Modified:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec
seismo/2D/SPECFEM2D/trunk/setup/constants.h.in
seismo/2D/SPECFEM2D/trunk/src/specfem2D/Makefile.in
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
adds Zhinan's pml routines
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec 2012-06-09 00:55:48 UTC (rev 20341)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec 2012-06-10 21:11:46 UTC (rev 20342)
@@ -73,7 +73,7 @@
USE_SNAPSHOT_NUMBER_IN_FILENAME = .false. # use snapshot number in the file name of JPEG color snapshots instead of the time step
#### for PostScript snapshots ####
-output_postscript_snapshot = .true. # output Postscript snapshot of the results
+output_postscript_snapshot = .false. # output Postscript snapshot of the results
imagetype_postscript = 1 # display 1=displ vector 2=veloc vector 3=accel vector; small arrows are displayed for the vectors
meshvect = .false. # display mesh on vector plots or not
modelvect = .false. # display velocity model on vector plots
Modified: seismo/2D/SPECFEM2D/trunk/setup/constants.h.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/setup/constants.h.in 2012-06-09 00:55:48 UTC (rev 20341)
+++ seismo/2D/SPECFEM2D/trunk/setup/constants.h.in 2012-06-10 21:11:46 UTC (rev 20342)
@@ -184,3 +184,7 @@
real(kind=CUSTOM_REAL), parameter :: r0_guenneau = 0.20_CUSTOM_REAL
real(kind=CUSTOM_REAL), parameter :: r1_guenneau = 0.40_CUSTOM_REAL
+! PML
+ logical, parameter :: PML_BOUNDARY_CONDITIONS = .false.
+ integer, parameter :: NELEM_PML_THICKNESS = 4
+
\ No newline at end of file
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/Makefile.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/Makefile.in 2012-06-09 00:55:48 UTC (rev 20341)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/Makefile.in 2012-06-10 21:11:46 UTC (rev 20342)
@@ -146,6 +146,7 @@
$O/paco_convolve_fft.o \
$O/plotgll.o \
$O/plotpost.o \
+ $O/pml_init.o \
$O/prepare_absorb.o \
$O/prepare_assemble_MPI.o \
$O/prepare_color_image.o \
@@ -357,6 +358,9 @@
$O/plotpost.o: ${S}/plotpost.F90 ${SETUP}/constants.h
${F90} $(FLAGS_CHECK) -c -o $O/plotpost.o ${S}/plotpost.F90
+$O/pml_init.o: ${S}/pml_init.F90 ${SETUP}/constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/pml_init.o ${S}/pml_init.F90
+
$O/prepare_absorb.o: ${S}/prepare_absorb.f90 ${SETUP}/constants.h
${F90} $(FLAGS_CHECK) -c -o $O/prepare_absorb.o ${S}/prepare_absorb.f90
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-06-09 00:55:48 UTC (rev 20341)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-06-10 21:11:46 UTC (rev 20342)
@@ -62,7 +62,12 @@
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k,&
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, &
- stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring)
+ stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring, &
+ is_PML,nspec_PML,npoin_PML,ibool_PML,spec_to_PML,which_PML_elem, &
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
+ rmemory_displ_elastic,rmemory_displ_elastic_corner, &
+ rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
+ rmemory_dux_dx_corner,rmemory_dux_dz_corner,rmemory_duz_dx_corner,rmemory_duz_dz_corner)
! compute forces for the elastic elements
@@ -185,8 +190,9 @@
real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl,jacobianl
! material properties of the elastic medium
- real(kind=CUSTOM_REAL) :: mul_unrelaxed_elastic,lambdal_unrelaxed_elastic,lambdaplus2mu_unrelaxed_elastic,kappal,cpl,csl,rhol, &
- lambdal_relaxed_viscoelastic,mul_relaxed_viscoelastic,lambdalplus2mul_relaxed_viscoel
+ real(kind=CUSTOM_REAL) :: mul_unrelaxed_elastic,lambdal_unrelaxed_elastic, &
+ lambdaplus2mu_unrelaxed_elastic,kappal,cpl,csl,rhol, &
+ lambdal_relaxed_viscoelastic,mul_relaxed_viscoelastic,lambdalplus2mul_relaxed_viscoel
! for attenuation
real(kind=CUSTOM_REAL) :: Un,Unp1,tauinv,Sn,Snp1,theta_n,theta_np1,tauinvsquare,tauinvcube,tauinvUn
@@ -209,6 +215,31 @@
integer :: ifirstelem,ilastelem
+!CPML coefficients and memory variables
+ integer :: nspec_PML,npoin_PML,iPML,ispec_PML
+ logical, dimension(4,nspec) :: which_PML_elem
+ logical, dimension(nspec) :: is_PML
+ integer, dimension(nspec) :: spec_to_PML
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool_PML
+
+ real(kind=CUSTOM_REAL), dimension(2,3,NGLLX,NGLLZ,nspec_PML) :: rmemory_displ_elastic,rmemory_displ_elastic_corner
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLZ,nspec_PML) :: &
+ rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz,&
+ rmemory_dux_dx_corner,rmemory_dux_dz_corner,rmemory_duz_dx_corner,rmemory_duz_dz_corner
+ real(kind=CUSTOM_REAL), dimension(npoin_PML) :: K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
+
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLZ,nspec_PML) :: accel_elastic_PML, accel_elastic_PML_corner
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) ::PML_dux_dxl,PML_dux_dzl,PML_duz_dxl,PML_duz_dzl,&
+ PML_dux_dxl_new,PML_dux_dzl_new,PML_duz_dxl_new,PML_duz_dzl_new
+ real(kind=CUSTOM_REAL) :: coef0_x, coef1_x, coef2_x, coef0_z, coef1_z, coef2_z,bb !,aa
+
+ real(kind=CUSTOM_REAL) :: A0, A1, A2, A3, A4, A5, A6, A7, A8 !, A9, A10
+
+ real(kind=CUSTOM_REAL) :: dux_dxi_new,dux_dgamma_new,duz_dxi_new,duz_dgamma_new !duy_dxi_new,duy_dgamma_new
+ real(kind=CUSTOM_REAL) :: dux_dxl_new,dux_dzl_new,duz_dxl_new,duz_dzl_new
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: displ_elastic_new
+
+
! this to avoid a warning at execution time about an undefined variable being used
! for the SH component in the case of a P-SV calculation, and vice versa
sigma_xx = 0
@@ -217,6 +248,21 @@
sigma_zy = 0
sigma_zz = 0
+ if( PML_BOUNDARY_CONDITIONS ) then
+ accel_elastic_PML = 0._CUSTOM_REAL
+ accel_elastic_PML_corner = 0._CUSTOM_REAL
+ PML_dux_dxl = 0._CUSTOM_REAL
+ PML_dux_dzl = 0._CUSTOM_REAL
+ PML_duz_dxl = 0._CUSTOM_REAL
+ PML_duz_dzl = 0._CUSTOM_REAL
+ PML_dux_dxl_new = 0._CUSTOM_REAL
+ PML_dux_dzl_new = 0._CUSTOM_REAL
+ PML_duz_dxl_new = 0._CUSTOM_REAL
+ PML_duz_dzl_new = 0._CUSTOM_REAL
+ displ_elastic_new = displ_elastic + deltat * veloc_elastic
+ endif
+
+
! compute Grad(displ_elastic) at time step n for attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
temp_displ_elastic = displ_elastic + deltat * veloc_elastic + 0.5 * deltat**2 * accel_elastic
@@ -323,6 +369,364 @@
duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
+ if( PML_BOUNDARY_CONDITIONS ) then
+ if(is_PML(ispec))then
+ ! derivative along x and along z
+ dux_dxi_new = ZERO
+ duz_dxi_new = ZERO
+ dux_dgamma_new = ZERO
+ duz_dgamma_new = ZERO
+
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the two loops because NGLLX == NGLLZ
+ do k = 1,NGLLX
+ dux_dxi_new = dux_dxi_new + displ_elastic_new(1,ibool(k,j,ispec))*hprime_xx(i,k)
+ duz_dxi_new = duz_dxi_new + displ_elastic_new(3,ibool(k,j,ispec))*hprime_xx(i,k)
+ dux_dgamma_new = dux_dgamma_new + displ_elastic_new(1,ibool(i,k,ispec))*hprime_zz(j,k)
+ duz_dgamma_new = duz_dgamma_new + displ_elastic_new(3,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)
+
+ ! derivatives of displacement
+ dux_dxl_new = dux_dxi_new*xixl + dux_dgamma_new*gammaxl
+ dux_dzl_new = dux_dxi_new*xizl + dux_dgamma_new*gammazl
+ duz_dxl_new = duz_dxi_new*xixl + duz_dgamma_new*gammaxl
+ duz_dzl_new = duz_dxi_new*xizl + duz_dgamma_new*gammazl
+ endif
+ endif
+
+ iglob = ibool(i,j,ispec)
+
+ if( PML_BOUNDARY_CONDITIONS ) then
+ if( is_PML(ispec) ) then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+
+!------------------------------------------------------------------------------
+!---------------------------- LEFT & RIGHT ------------------------------------
+!------------------------------------------------------------------------------
+ if ( which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec) ) then
+
+ !---------------------- A8 and A6 --------------------------
+ A8 = - d_x_store(iPML) / ( k_x_store(iPML) ** 2 )
+ A6 = 0.d0
+
+ bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
+ coef0_x = exp(-bb * deltat)
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = (1.d0 - exp(-bb * deltat / 2.d0)) * A8 / bb
+ coef2_x = (1.d0 - exp(-bb* deltat / 2.d0)) * exp(-bb * deltat / 2.d0) * A8 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A8
+ coef2_x = deltat * A8
+ end if
+
+ PML_dux_dxl(i,j,ispec_PML) = dux_dxl
+ PML_dux_dxl_new(i,j,ispec_PML) = dux_dxl_new
+
+ rmemory_dux_dx(1,i,j,ispec_PML) = 0.d0
+ !! DK DK new from Wang eq (21)
+ rmemory_dux_dx(2,i,j,ispec_PML) = coef0_x*rmemory_dux_dx(2,i,j,ispec_PML) &
+ + PML_dux_dxl_new(i,j,ispec_PML) * coef1_x + PML_dux_dxl(i,j,ispec_PML) * coef2_x
+
+ PML_duz_dxl(i,j,ispec_PML)=duz_dxl !!!when in corner this is no need
+ PML_duz_dxl_new(i,j,ispec_PML)=duz_dxl_new !!!when in corner this is no need
+
+ rmemory_duz_dx(1,i,j,ispec_PML) = 0.d0
+ !! DK DK new from Wang eq (21)
+ rmemory_duz_dx(2,i,j,ispec_PML) = coef0_x*rmemory_duz_dx(2,i,j,ispec_PML) &
+ + PML_duz_dxl_new(i,j,ispec_PML) * coef1_x + PML_duz_dxl(i,j,ispec_PML) * coef2_x
+
+ !---------------------- A5 and A7 --------------------------
+ A5 = d_x_store(iPML )
+ A7 = 0.d0
+
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs( bb ) > 1.d-3) then
+ coef1_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * A5 / bb
+ coef2_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) * &
+ A5 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A5
+ coef2_x = deltat * A5 ! Rene Matzen
+ end if
+
+ PML_dux_dzl(i,j,ispec_PML)=dux_dzl !!!when in corner this is no need
+ PML_dux_dzl_new(i,j,ispec_PML)=dux_dzl_new !!!when in corner this is no need
+
+ !! DK DK new from Wang eq (21)
+ rmemory_dux_dz(1,i,j,ispec_PML) = coef0_x * rmemory_dux_dz(1,i,j,ispec_PML) &
+ + PML_dux_dzl_new(i,j,ispec_PML) *coef1_x + PML_dux_dzl(i,j,ispec_PML) * coef2_x
+ rmemory_dux_dz(2,i,j,ispec_PML) = 0.d0
+
+ PML_duz_dzl(i,j,ispec_PML)=duz_dzl !!!when in corner this is no need
+ PML_duz_dzl_new(i,j,ispec_PML)=duz_dzl_new !!!when in corner this is no need
+
+ !! DK DK new from Wang eq (21)
+ rmemory_duz_dz(1,i,j,ispec_PML) = coef0_x * rmemory_duz_dz(1,i,j,ispec_PML) &
+ + PML_duz_dzl_new(i,j,ispec_PML) *coef1_x + PML_duz_dzl(i,j,ispec_PML) * coef2_x
+ rmemory_duz_dz(2,i,j,ispec_PML) = 0.d0
+
+ dux_dxl = PML_dux_dxl(i,j,ispec_PML) + rmemory_dux_dx(1,i,j,ispec_PML) + rmemory_dux_dx(2,i,j,ispec_PML)
+ duz_dxl = PML_duz_dxl(i,j,ispec_PML) + rmemory_duz_dx(1,i,j,ispec_PML) + rmemory_duz_dx(2,i,j,ispec_PML)
+ dux_dzl = PML_dux_dzl(i,j,ispec_PML) + rmemory_dux_dz(1,i,j,ispec_PML) + rmemory_dux_dz(2,i,j,ispec_PML)
+ duz_dzl = PML_duz_dzl(i,j,ispec_PML) + rmemory_duz_dz(1,i,j,ispec_PML) + rmemory_duz_dz(2,i,j,ispec_PML)
+
+!------------------------------------------------------------------------------
+!---------------------------- CORNER ------------------------------------------
+!------------------------------------------------------------------------------
+ if ( which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec) ) then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+
+ !----------------------- A7 ------------------------------
+ if ( abs( alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A7 = 0.d0
+ elseif ( abs( alpha_x_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) ) < 1.d-3 &
+ ) then
+ A7 = d_z_store(iPML)
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) - d_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( k_x_store(iPML) - k_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A7 = 0.d0
+ else
+ A7 = d_z_store(iPML) * ( alpha_x_store(iPML) - alpha_z_store(iPML) ) &
+ / ( k_x_store(iPML) * ( alpha_x_store(iPML) - alpha_z_store(iPML) ) + d_x_store(iPML) )
+ end if
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = ( 1 - exp(- bb * deltat / 2) ) * A7 / bb
+ coef2_x = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) * A7 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A7
+ coef2_x = deltat * A7 ! Rene Matzen
+ end if
+
+ !! DK DK new from Wang eq (21)
+ rmemory_dux_dx_corner(1,i,j,ispec_PML) = coef0_x*rmemory_dux_dx_corner(1,i,j,ispec_PML) &
+ + PML_dux_dxl_new(i,j,ispec_PML) * coef1_x + PML_dux_dxl(i,j,ispec_PML) * coef2_x
+
+ !! DK DK new from Wang eq (21)
+ rmemory_duz_dx_corner(1,i,j,ispec_PML) = coef0_x*rmemory_duz_dx_corner(1,i,j,ispec_PML) &
+ + PML_duz_dxl_new(i,j,ispec_PML) * coef1_x + PML_duz_dxl(i,j,ispec_PML) * coef2_x
+
+ !---------------------------- A8 ----------------------------
+ if ( abs( alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A8 = - d_x_store(iPML) / ( k_x_store(iPML) ** 2 )
+ elseif ( abs( alpha_x_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) ) < 1.d-3 &
+ ) then
+ A8 = 0.d0
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) - d_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( k_x_store(iPML) - k_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A8 = 0.d0
+ else
+ A8 = d_x_store(iPML) * &
+ ( k_x_store(iPML) * k_z_store(iPML) * ( alpha_z_store(iPML) - alpha_x_store(iPML) ) &
+ + k_x_store(iPML) * d_z_store(iPML) - k_z_store(iPML) * d_x_store(iPML) &
+ ) &
+ / ( k_x_store(iPML)**2 * ( k_x_store(iPML) * ( alpha_x_store(iPML) - alpha_z_store(iPML) ) &
+ + d_x_store(iPML) ) &
+ )
+ end if
+ bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
+ coef0_z = exp(- bb * deltat)
+
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_z = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * A8 / bb
+ coef2_z = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * exp(- bb * deltat / 2.d0) * A8 / bb
+ else
+ coef1_z = deltat / 2.0d0 * A8
+ coef2_z = deltat * A8 ! Rene Matzen
+ end if
+
+ !! DK DK new from Wang eq (21)
+ rmemory_dux_dx_corner(2,i,j,ispec_PML) = coef0_z*rmemory_dux_dx_corner(2,i,j,ispec_PML) &
+ + PML_dux_dxl_new(i,j,ispec_PML) * coef1_z + PML_dux_dxl(i,j,ispec_PML) * coef2_z
+
+ !! DK DK new from Wang eq (21)
+ rmemory_duz_dx_corner(2,i,j,ispec_PML) = coef0_z*rmemory_duz_dx_corner(2,i,j,ispec_PML) &
+ + PML_duz_dxl_new(i,j,ispec_PML) * coef1_z + PML_duz_dxl(i,j,ispec_PML) * coef2_z
+
+ !---------------------------- A5 ----------------------------
+ if ( abs( alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A5 = d_x_store(iPML)
+ elseif ( abs( alpha_x_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) ) < 1.d-3 &
+ ) then
+ A5 = 0.d0
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) - d_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( k_x_store(iPML) - k_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A5 = 0.d0
+ else
+ A5 = d_x_store(iPML) * ( alpha_z_store(iPML) - alpha_x_store(iPML) ) &
+ / ( k_z_store(iPML) * ( alpha_z_store(iPML) - alpha_x_store(iPML) ) + d_z_store(iPML) )
+ end if
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = ( 1 - exp(- bb * deltat / 2) ) * A5 / bb
+ coef2_x = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) * A5 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A5
+ coef2_x = deltat * A5 ! Rene Matzen
+ end if
+
+ !! DK DK new from Wang eq (21)
+ rmemory_dux_dz_corner(1,i,j,ispec_PML) = coef0_x * rmemory_dux_dz_corner(1,i,j,ispec_PML) &
+ + PML_dux_dzl_new(i,j,ispec_PML) *coef1_x + PML_dux_dzl(i,j,ispec_PML) * coef2_x
+
+ !! DK DK new from Wang eq (21)
+ rmemory_duz_dz_corner(1,i,j,ispec_PML) = coef0_x * rmemory_duz_dz_corner(1,i,j,ispec_PML) &
+ + PML_duz_dzl_new(i,j,ispec_PML) *coef1_x + PML_duz_dzl(i,j,ispec_PML) * coef2_x
+
+ !---------------------------- A6 ----------------------------
+ if ( abs( alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A6 = 0.d0
+ elseif ( abs( alpha_x_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) ) < 1.d-3 &
+ ) then
+ A6 = - d_z_store(iPML) / ( k_z_store(iPML) ** 2 )
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) - d_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( k_x_store(iPML) - k_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A6 = 0.d0
+ else
+ A6 = d_z_store(iPML) * &
+ ( k_z_store(iPML) * k_x_store(iPML) * ( alpha_x_store(iPML) - alpha_z_store(iPML) ) &
+ + k_z_store(iPML) * d_x_store(iPML) - k_x_store(iPML) * d_z_store(iPML) &
+ ) &
+ / ( k_z_store(iPML)**2 * ( k_z_store(iPML) * ( alpha_z_store(iPML) - alpha_x_store(iPML) ) &
+ + d_z_store(iPML) ) &
+ )
+ end if
+ bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
+ coef0_z = exp(- bb * deltat)
+
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_z = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * A6 / bb
+ coef2_z = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * exp(- bb * deltat / 2.d0) * A6 / bb
+ else
+ coef1_z = deltat / 2.0d0 * A6
+ coef2_z = deltat * A6 ! Rene Matzen
+ end if
+
+ !! DK DK new from Wang eq (21)
+ rmemory_dux_dz_corner(2,i,j,ispec_PML) = coef0_z * rmemory_dux_dz_corner(2,i,j,ispec_PML) &
+ + PML_dux_dzl_new(i,j,ispec_PML) *coef1_z + PML_dux_dzl(i,j,ispec_PML) * coef2_z
+
+ !! DK DK new from Wang eq (21)
+ rmemory_duz_dz_corner(2,i,j,ispec_PML) = coef0_z * rmemory_duz_dz_corner(2,i,j,ispec_PML) &
+ + PML_duz_dzl_new(i,j,ispec_PML) *coef1_z + PML_duz_dzl(i,j,ispec_PML) * coef2_z
+
+ dux_dxl = PML_dux_dxl(i,j,ispec_PML) + rmemory_dux_dx_corner(1,i,j,ispec_PML) +&
+ rmemory_dux_dx_corner(2,i,j,ispec_PML)
+ duz_dxl = PML_duz_dxl(i,j,ispec_PML) + rmemory_duz_dx_corner(1,i,j,ispec_PML) +&
+ rmemory_duz_dx_corner(2,i,j,ispec_PML)
+ dux_dzl = PML_dux_dzl(i,j,ispec_PML) + rmemory_dux_dz_corner(1,i,j,ispec_PML) + &
+ rmemory_dux_dz_corner(2,i,j,ispec_PML)
+ duz_dzl = PML_duz_dzl(i,j,ispec_PML) + rmemory_duz_dz_corner(1,i,j,ispec_PML) + &
+ rmemory_duz_dz_corner(2,i,j,ispec_PML)
+ endif
+
+ else
+!------------------------------------------------------------------------------
+!---------------------------- TOP & BOTTOM ------------------------------------
+!------------------------------------------------------------------------------
+
+ !---------------------- A5 and A7 --------------------------
+ A7 = d_z_store(iPML )
+ A5 = 0.d0
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs( bb ) > 1.d-3) then
+ coef1_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * A7 / bb
+ coef2_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) * &
+ A7 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A7
+ coef2_x = deltat * A7 ! Rene Matzen
+ end if
+
+ PML_dux_dxl(i,j,ispec_PML) = dux_dxl !!!when in corner this is no need
+ PML_dux_dxl_new(i,j,ispec_PML) = dux_dxl_new !!!when in corner this is no need
+
+ !! DK DK new from Wang eq (21)
+ rmemory_dux_dx(1,i,j,ispec_PML) = coef0_x*rmemory_dux_dx(1,i,j,ispec_PML) &
+ + PML_dux_dxl_new(i,j,ispec_PML) * coef1_x + PML_dux_dxl(i,j,ispec_PML) * coef2_x
+ rmemory_dux_dx(2,i,j,ispec_PML) = 0.d0
+
+ PML_duz_dxl(i,j,ispec_PML)=duz_dxl !!!when in corner this is no need
+ PML_duz_dxl_new(i,j,ispec_PML)=duz_dxl_new !!!when in corner this is no need
+
+ !! DK DK new from Wang eq (21)
+ rmemory_duz_dx(1,i,j,ispec_PML) = coef0_x*rmemory_duz_dx(1,i,j,ispec_PML) &
+ + PML_duz_dxl_new(i,j,ispec_PML) * coef1_x + PML_duz_dxl(i,j,ispec_PML) * coef2_x
+ rmemory_duz_dx(2,i,j,ispec_PML) = 0.d0
+
+ !---------------------- A8 and A6 --------------------------
+ A6 = - d_z_store(iPML) / ( k_z_store(iPML) ** 2 )
+ A8 = 0.d0
+ bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
+ coef0_x = exp(-bb * deltat)
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = (1.d0 - exp(-bb * deltat / 2.d0)) * A6 / bb
+ coef2_x = (1.d0 - exp(-bb* deltat / 2.d0)) * exp(-bb * deltat / 2.d0) * A6 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A6
+ coef2_x = deltat * A6
+ end if
+
+ PML_dux_dzl(i,j,ispec_PML)=dux_dzl !!!when in corner this is no need
+ PML_dux_dzl_new(i,j,ispec_PML)=dux_dzl_new !!!when in corner this is no need
+
+ rmemory_dux_dz(1,i,j,ispec_PML) = 0.d0
+ !! DK DK new from Wang eq (21)
+ rmemory_dux_dz(2,i,j,ispec_PML) = coef0_x * rmemory_dux_dz(2,i,j,ispec_PML) &
+ + PML_dux_dzl_new(i,j,ispec_PML) *coef1_x + PML_dux_dzl(i,j,ispec_PML) * coef2_x
+
+ PML_duz_dzl(i,j,ispec_PML)=duz_dzl !!!when in corner this is no need
+ PML_duz_dzl_new(i,j,ispec_PML)=duz_dzl_new !!!when in corner this is no need
+
+ rmemory_duz_dz(1,i,j,ispec_PML) = 0.d0
+ !! DK DK new from Wang eq (21)
+ rmemory_duz_dz(2,i,j,ispec_PML) = coef0_x * rmemory_duz_dz(2,i,j,ispec_PML) &
+ + PML_duz_dzl_new(i,j,ispec_PML) *coef1_x + PML_duz_dzl(i,j,ispec_PML) * coef2_x
+
+ dux_dxl = dux_dxl + rmemory_dux_dx(1,i,j,ispec_PML) + rmemory_dux_dx(2,i,j,ispec_PML)
+ duz_dxl = duz_dxl + rmemory_duz_dx(1,i,j,ispec_PML) + rmemory_duz_dx(2,i,j,ispec_PML)
+ dux_dzl = dux_dzl + rmemory_dux_dz(1,i,j,ispec_PML) + rmemory_dux_dz(2,i,j,ispec_PML)
+ duz_dzl = duz_dzl + rmemory_duz_dz(1,i,j,ispec_PML) + rmemory_duz_dz(2,i,j,ispec_PML)
+
+ endif
+ endif
+ endif ! PML_BOUNDARY_CONDITIONS
+
if(SIMULATION_TYPE == 2) then ! Adjoint calculation, backward wavefield
b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
@@ -390,6 +794,7 @@
sigma_xz = mul_unrelaxed_elastic*(duz_dxl + dux_dzl)
sigma_zy = mul_unrelaxed_elastic*duy_dzl
sigma_zz = lambdaplus2mu_unrelaxed_elastic*duz_dzl + lambdal_unrelaxed_elastic*dux_dxl
+ sigma_zx = sigma_xz
!! DK DK added this for Guenneau, March 2012
#ifdef USE_GUENNEAU
@@ -397,6 +802,17 @@
#endif
!! DK DK added this for Guenneau, March 2012
+ if( PML_BOUNDARY_CONDITIONS ) then
+ if(is_PML(ispec)) then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ sigma_xx = lambdaplus2mu_unrelaxed_elastic*dux_dxl + lambdal_unrelaxed_elastic*PML_duz_dzl(i,j,ispec_PML)
+ sigma_zz = lambdaplus2mu_unrelaxed_elastic*duz_dzl + lambdal_unrelaxed_elastic*PML_dux_dxl(i,j,ispec_PML)
+ sigma_zx = mul_unrelaxed_elastic * (PML_duz_dxl(i,j,ispec_PML) + dux_dzl)
+ sigma_xz = mul_unrelaxed_elastic * (PML_dux_dzl(i,j,ispec_PML) + duz_dxl)
+ endif
+ endif
+
if(SIMULATION_TYPE == 2) then ! Adjoint calculation, backward wavefield
b_sigma_xx = lambdaplus2mu_unrelaxed_elastic*b_dux_dxl + lambdal_unrelaxed_elastic*b_duz_dzl
b_sigma_xy = mul_unrelaxed_elastic*b_duy_dxl
@@ -458,7 +874,7 @@
jacobianl = jacobian(i,j,ispec)
! the stress tensor is symmetric by default, unless defined otherwise
-#ifndef USE_GUENNEAU
+#ifdef USE_GUENNEAU
sigma_zx = sigma_xz
#endif
@@ -485,6 +901,243 @@
enddo
enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! PML
+ if( PML_BOUNDARY_CONDITIONS ) then
+ if(is_PML(ispec)) then
+ ! first double loop over GLL points to compute and store gradients
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ if ( assign_external_model) then
+ rhol = rhoext(i,j,ispec)
+ else
+ rhol = density(1,kmato(ispec))
+ endif
+
+ if ( is_PML( ispec ) ) then
+ iPML=ibool_PML(i,j,ispec)
+ iglob=ibool(i,j,ispec)
+
+!------------------------------------------------------------------------------
+!---------------------------- LEFT & RIGHT ------------------------------------
+!------------------------------------------------------------------------------
+ if ( which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec) ) then
+
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ iglob=ibool(i,j,ispec)
+
+ !---------------------- A3 and A4 --------------------------
+ A3 = d_x_store(iPML ) * alpha_x_store(iPML) ** 2
+ A4 = 0.d0
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs( bb ) > 1.d-3) then
+ coef1_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * A3 / bb
+ coef2_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) * &
+ A3 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A3
+ coef2_x = deltat * A3 ! Rene Matzen
+ end if
+
+ rmemory_displ_elastic(1,1,i,j,ispec_PML)=coef0_x * rmemory_displ_elastic(1,1,i,j,ispec_PML) &
+ + displ_elastic_new(1,iglob) * coef1_x + displ_elastic(1,iglob) * coef2_x
+ rmemory_displ_elastic(2,1,i,j,ispec_PML)=0.0
+
+ rmemory_displ_elastic(1,3,i,j,ispec_PML)=coef0_x * rmemory_displ_elastic(1,3,i,j,ispec_PML) &
+ + displ_elastic_new(3,iglob) * coef1_x + displ_elastic(3,iglob) * coef2_x
+ rmemory_displ_elastic(2,3,i,j,ispec_PML)=0.0
+
+!------------------------------------------------------------------------------
+!-------------------------------- CORNER --------------------------------------
+!------------------------------------------------------------------------------
+ if ( which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec) ) then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+
+ !---------------------------- A3 ----------------------------
+ if ( abs( alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A3 = d_x_store(iPML) * alpha_x_store(iPML) ** 2
+ elseif ( abs( alpha_x_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) ) < 1.d-3 &
+ ) then
+ A3 = 0.d0
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) - d_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( k_x_store(iPML) - k_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A3 = 0.d0
+ else
+ A3 = alpha_x_store(iPML) ** 2 * d_x_store(iPML) * &
+ ( k_z_store(iPML) * ( alpha_x_store(iPML) - alpha_z_store(iPML) ) - d_z_store(iPML) ) &
+ / ( alpha_x_store(iPML) - alpha_z_store(iPML) )
+ end if
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = ( 1 - exp(- bb * deltat / 2) ) * A3 / bb
+ coef2_x = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) * A3 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A3
+ coef2_x = deltat * A3 ! Rene Matzen
+ end if
+
+ rmemory_displ_elastic_corner(1,1,i,j,ispec_PML)= &
+ coef0_x * rmemory_displ_elastic_corner(1,1,i,j,ispec_PML) &
+ + displ_elastic_new(1,iglob) * coef1_x + displ_elastic(1,iglob) * coef2_x
+ rmemory_displ_elastic_corner(1,3,i,j,ispec_PML)= &
+ coef0_x * rmemory_displ_elastic_corner(1,3,i,j,ispec_PML) &
+ + displ_elastic_new(3,iglob) * coef1_x + displ_elastic(3,iglob) * coef2_x
+
+
+ !---------------------------- A4 ----------------------------
+ if ( abs( alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A4 = 0.d0
+ elseif ( abs( alpha_x_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) ) < 1.d-3 &
+ ) then
+ A4 = d_z_store(iPML) * alpha_z_store(iPML) ** 2
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) - d_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( k_x_store(iPML) - k_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A4 = 0.d0
+ else
+ A4 = alpha_z_store(iPML) ** 2 * d_z_store(iPML) * &
+ ( k_x_store(iPML) * ( alpha_z_store(iPML) - alpha_x_store(iPML) ) - d_x_store(iPML) ) &
+ / ( alpha_z_store(iPML) - alpha_x_store(iPML) )
+ end if
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = ( 1 - exp(- bb * deltat / 2) ) * A4 / bb
+ coef2_x = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) * A4 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A4
+ coef2_x = deltat * A4 ! Rene Matzen
+ end if
+
+ rmemory_displ_elastic_corner(2,1,i,j,ispec_PML)= &
+ coef0_x * rmemory_displ_elastic_corner(2,1,i,j,ispec_PML) &
+ + displ_elastic_new(1,iglob) * coef1_x + displ_elastic(1,iglob) * coef2_x
+ rmemory_displ_elastic_corner(2,3,i,j,ispec_PML)= &
+ coef0_x * rmemory_displ_elastic_corner(2,3,i,j,ispec_PML) &
+ + displ_elastic_new(3,iglob) * coef1_x + displ_elastic(3,iglob) * coef2_x
+
+ endif
+
+ else
+!------------------------------------------------------------------------------
+!-------------------------------- TOP & BOTTOM --------------------------------
+!------------------------------------------------------------------------------
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ iglob=ibool(i,j,ispec)
+
+ !---------------------- A3 and A4 ----------------------------
+ A4 = d_z_store(iPML ) * alpha_z_store(iPML) ** 2
+ A3 = 0.d0
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs( bb ) > 1.d-3) then
+ coef1_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * A4 / bb
+ coef2_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) * &
+ A4 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A4
+ coef2_x = deltat * A4 ! Rene Matzen
+ end if
+
+ rmemory_displ_elastic(1,1,i,j,ispec_PML)=0.d0
+ rmemory_displ_elastic(2,1,i,j,ispec_PML)=coef0_x * rmemory_displ_elastic(2,1,i,j,ispec_PML) &
+ + displ_elastic_new(1,iglob) * coef1_x + displ_elastic(1,iglob) * coef2_x
+
+ rmemory_displ_elastic(1,3,i,j,ispec_PML)=0.d0
+ rmemory_displ_elastic(2,3,i,j,ispec_PML)=coef0_x * rmemory_displ_elastic(2,3,i,j,ispec_PML) &
+ + displ_elastic_new(3,iglob) * coef1_x + displ_elastic(3,iglob) * coef2_x
+
+ endif
+
+
+ if ( which_PML_elem(ILEFT,ispec) .or. which_PML_elem(IRIGHT,ispec)) then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ iglob=ibool(i,j,ispec)
+
+ A0 = - alpha_x_store( iPML ) * d_x_store(iPML) * k_z_store(iPML)
+ A1 = d_x_store(iPML) * k_z_store(iPML)
+ A2 = k_x_store(iPML) * k_z_store(iPML)
+
+ accel_elastic_PML(1,i,j,ispec_PML) = wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
+ ( A1 * veloc_elastic(1,iglob)+ &
+ rmemory_displ_elastic(1,1,i,j,ispec_PML)+rmemory_displ_elastic(2,1,i,j,ispec_PML)&
+ + A0 * displ_elastic(1,iglob) )
+ accel_elastic_PML(3,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
+ ( A1 * veloc_elastic(3,iglob)+&
+ rmemory_displ_elastic(1,3,i,j,ispec_PML)+rmemory_displ_elastic(2,3,i,j,ispec_PML)&
+ + A0 *displ_elastic(3,iglob) )
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
+ if ( which_PML_elem(ITOP,ispec) .or. which_PML_elem(IBOTTOM,ispec)) then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+
+ A0 = d_x_store(iPML) * d_z_store(iPML) &
+ - alpha_x_store( iPML ) * d_x_store(iPML) * k_z_store(iPML) &
+ - alpha_z_store( iPML ) * d_z_store(iPML) * k_x_store(iPML)
+
+ A1 = d_x_store(iPML) * k_z_store(iPML) + d_z_store(iPML) * k_x_store(iPML)
+
+ A2 = k_x_store(iPML) * k_z_store(iPML)
+
+ accel_elastic_PML_corner(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
+ ( A1 *veloc_elastic(1,iglob)+ &
+ rmemory_displ_elastic_corner(1,1,i,j,ispec_PML)+rmemory_displ_elastic_corner(2,1,i,j,ispec_PML)&
+ + A0 * displ_elastic(1,iglob) )
+ accel_elastic_PML_corner(3,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
+ ( A1 * veloc_elastic(3,iglob)+ &
+ rmemory_displ_elastic_corner(1,3,i,j,ispec_PML)+rmemory_displ_elastic_corner(2,3,i,j,ispec_PML)&
+ + A0 * displ_elastic(3,iglob) )
+
+ endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
+ else
+
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ iglob=ibool(i,j,ispec)
+
+ A0 = - alpha_z_store( iPML ) * d_z_store(iPML) * k_x_store(iPML)
+ A1 = d_z_store(iPML) * k_x_store(iPML)
+ A2 = k_x_store(iPML) * k_z_store(iPML)
+
+ accel_elastic_PML(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
+ ( A1 * veloc_elastic(1,iglob)+ &
+ rmemory_displ_elastic(1,1,i,j,ispec_PML)+rmemory_displ_elastic(2,1,i,j,ispec_PML)&
+ + A0 * displ_elastic(1,iglob) )
+ accel_elastic_PML(3,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
+ ( A1 * veloc_elastic(3,iglob)+&
+ rmemory_displ_elastic(1,3,i,j,ispec_PML)+rmemory_displ_elastic(2,3,i,j,ispec_PML)&
+ + A0 * displ_elastic(3,iglob) )
+ endif
+ endif
+
+ enddo
+ enddo
+
+ endif ! end of test if elastic element
+ endif ! PML_BOUNDARY_CONDITIONS
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
!
! second double-loop over GLL to compute all the terms
!
@@ -497,9 +1150,12 @@
! and assemble the contributions
! 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) - (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))
+ 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) &
+ - (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(SIMULATION_TYPE == 2) then ! Adjoint calculation, backward wavefield
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
@@ -511,6 +1167,42 @@
endif
enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if( PML_BOUNDARY_CONDITIONS ) then
+ if(is_PML(ispec))then
+ if (which_PML_elem(ILEFT,ispec) .or. which_PML_elem(IRIGHT,ispec)) then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - accel_elastic_PML(1,i,j,ispec_PML)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - accel_elastic_PML(3,i,j,ispec_PML)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
+ if (which_PML_elem(ITOP,ispec) .or. which_PML_elem(IBOTTOM,ispec)) then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + accel_elastic_PML(1,i,j,ispec_PML)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + accel_elastic_PML(3,i,j,ispec_PML)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - accel_elastic_PML_corner(1,i,j,ispec_PML)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - accel_elastic_PML_corner(3,i,j,ispec_PML)
+ endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
+ else
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - accel_elastic_PML(1,i,j,ispec_PML)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - accel_elastic_PML(3,i,j,ispec_PML)
+ endif
+ if(it==900)then
+ write(IOUT,*)accel_elastic_PML(1,i,j,ispec_PML),accel_elastic_PML(3,i,j,ispec_PML),'accel_elastic_PML'
+ endif
+ endif
+ endif ! PML_BOUNDARY_CONDITIONS
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
enddo ! second loop over the GLL points
enddo
@@ -518,6 +1210,8 @@
enddo ! end of loop over all spectral elements
+
+
!
!--- absorbing boundaries
!
@@ -1023,10 +1717,13 @@
endif
elseif(SIMULATION_TYPE == 2) then
if(p_sv)then !P-SV waves
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_top(1,i,ib_top(ispecabs),NSTEP-it+1)
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_top(3,i,ib_top(ispecabs),NSTEP-it+1)
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) &
+ - b_absorb_elastic_top(1,i,ib_top(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) &
+ - b_absorb_elastic_top(3,i,ib_top(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_top(ispecabs),NSTEP-it+1)
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) &
+ - b_absorb_elastic_top(2,i,ib_top(ispecabs),NSTEP-it+1)
endif
endif
@@ -1040,6 +1737,72 @@
endif ! end of absorbing boundaries
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ !--- absorbing boundaries
+ !
+ if( PML_BOUNDARY_CONDITIONS ) then
+
+! we have to put Dirichlet on the boundary of the PML
+ do ispecabs = 1,nelemabs
+
+ ispec = numabs(ispecabs)
+
+ if (is_PML(ispec)) then
+ ispec_PML=spec_to_PML(ispec)
+!--- left absorbing boundary
+ if(codeabs(ILEFT,ispecabs)) then
+ i = 1
+ do j = 1,NGLLZ
+ iglob = ibool(i,j,ispec)
+ displ_elastic(:,iglob) = 0._CUSTOM_REAL
+ veloc_elastic(:,iglob) = 0._CUSTOM_REAL
+ accel_elastic(:,iglob) = 0._CUSTOM_REAL
+ enddo
+ endif
+!--- right absorbing boundary
+ if(codeabs(IRIGHT,ispecabs)) then
+ i = NGLLX
+ do j = 1,NGLLZ
+ iglob = ibool(i,j,ispec)
+ displ_elastic(:,iglob) = 0._CUSTOM_REAL
+ veloc_elastic(:,iglob) = 0._CUSTOM_REAL
+ accel_elastic(:,iglob) = 0._CUSTOM_REAL
+ enddo
+
+ endif
+!--- bottom absorbing boundary
+ if(codeabs(IBOTTOM,ispecabs)) then
+ j = 1
+! exclude corners to make sure there is no contradiction on the normal
+ ibegin = 1
+ iend = NGLLX
+ do i = ibegin,iend
+ iglob = ibool(i,j,ispec)
+ displ_elastic(:,iglob) = 0._CUSTOM_REAL
+ veloc_elastic(:,iglob) = 0._CUSTOM_REAL
+ accel_elastic(:,iglob) = 0._CUSTOM_REAL
+ enddo
+ endif
+!--- top absorbing boundary
+ if(codeabs(ITOP,ispecabs)) then
+ j = NGLLZ
+! exclude corners to make sure there is no contradiction on the normal
+ ibegin = 1
+ iend = NGLLX
+ do i = ibegin,iend
+ iglob = ibool(i,j,ispec)
+ displ_elastic(:,iglob) = 0._CUSTOM_REAL
+ veloc_elastic(:,iglob) = 0._CUSTOM_REAL
+ accel_elastic(:,iglob) = 0._CUSTOM_REAL
+ enddo
+ endif ! end of top absorbing boundary
+ endif ! end of is_PML
+ enddo ! end specabs loop
+ endif ! end of absorbing boundaries PML_BOUNDARY_CONDITIONS
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
! --- add the source if it is a moment tensor
if(.not. initialfield) then
@@ -1172,7 +1935,8 @@
e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
endif
- e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + weight_rk * e1_force_RK(i,j,ispec,i_sls,i_stage)
+ e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) &
+ + weight_rk * e1_force_RK(i,j,ispec,i_sls,i_stage)
elseif(i_stage==4)then
@@ -1222,7 +1986,8 @@
if(i_stage==1)then
e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
endif
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) &
+ + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
elseif(i_stage==4)then
e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
(e11_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e11_force_RK(i,j,ispec,i_sls,2) + &
@@ -1267,7 +2032,8 @@
if(i_stage==1)then
e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
endif
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) &
+ + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
elseif(i_stage==4)then
e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
(e13_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e13_force_RK(i,j,ispec,i_sls,2) + &
Added: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-06-10 21:11:46 UTC (rev 20342)
@@ -0,0 +1,613 @@
+!========================================================================
+!
+! S P E C F E M 2 D Version 6 . 2
+! ------------------------------
+!
+! Copyright Universite de Pau, CNRS and INRIA, France,
+! and Princeton University / California Institute of Technology, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+ subroutine pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
+ nspec_PML,is_PML,which_PML_elem,which_PML_poin,spec_to_PML,ibool_PML, &
+ npoin_PML,icorner_iglob)
+
+ implicit none
+ include 'constants.h'
+
+#ifdef USE_MPI
+ include 'mpif.h'
+#endif
+
+ integer :: nspec,nglob,nelemabs,nspec_PML,npoin_PML
+ logical :: anyabs
+
+ integer :: ibound,ispecabs,ncorner,ispec,iglob
+ integer :: i,j,k,i_coef
+
+ logical, dimension(4,nspec) :: which_PML_elem
+ logical, dimension(4,nglob) :: which_PML_poin
+ integer, dimension(nglob) :: icorner_iglob
+ integer, dimension(nelemabs) :: numabs
+ logical, dimension(4,nelemabs) :: codeabs
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+ integer, dimension(nspec) :: spec_to_PML
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool_PML
+ integer, dimension(:), allocatable :: iPML_to_iglob
+ logical, dimension(nspec) :: is_PML
+
+ !!!!RM RM RM Beginning PML implementation of coefficients properties
+ ! allocate (is_PML(nspec))
+ ! is_PML(:)=.false.
+
+ nspec_PML=0
+
+ !!!detection of PML elements
+
+ !ibound is the side we are looking (bottom, right, top or left)
+ do ibound=1,4
+
+ icorner_iglob = ZERO
+ ncorner=0
+
+ if (anyabs) then
+ !mark any elements on the boundary as PML and list their corners
+ do ispecabs = 1,nelemabs
+ ispec = numabs(ispecabs)
+
+ !array to know which PML it is
+ which_PML_elem(ibound,ispec)=codeabs(ibound,ispecabs)
+
+ if(codeabs(ibound,ispecabs)) then ! we are on the good absorbing boundary
+ do j=1,NGLLZ,NGLLZ-1
+ do i=1,NGLLX,NGLLX-1
+ iglob=ibool(i,j,ispec)
+ k=1
+ do while(k<=ncorner .and. icorner_iglob(k)/=iglob)
+ k=k+1
+ end do
+ ncorner=ncorner+1
+ icorner_iglob(ncorner) = iglob
+
+ !array to know which PML it is
+ which_PML_poin(ibound,iglob) = codeabs(ibound,ispecabs)
+
+ enddo
+ enddo
+ endif ! we are on the good absorbing boundary
+
+ enddo
+ end if
+
+ !find elements stuck to boundary elements to define the 4 elements PML thickness
+ !we take 4 elements for the PML thickness
+ do i_coef=2,NELEM_PML_THICKNESS
+
+ do ispec=1,nspec
+ if (.not.which_PML_elem(ibound,ispec)) then
+ do j=1,NGLLZ,NGLLZ-1
+ do i=1,NGLLX,NGLLX-1
+ iglob=ibool(i,j,ispec)
+ do k=1,ncorner
+ if (iglob==icorner_iglob(k)) then
+ which_PML_elem(ibound,ispec) = .true.
+ end if
+ end do
+ end do
+ end do
+ end if
+ end do
+
+ !list every corner of each PML elements detected
+ ncorner=0
+ icorner_iglob=ZERO
+ nspec_PML=0
+ do ispec=1,nspec
+ if (which_PML_elem(ibound,ispec)) then
+ is_PML(ispec)=.true.
+ do j=1,NGLLZ,NGLLZ-1
+ do i=1,NGLLX,NGLLX-1
+ iglob=ibool(i,j,ispec)
+ k=1
+ do while(k<=ncorner .and. icorner_iglob(k)/=iglob)
+ k=k+1
+ end do
+ ncorner=ncorner+1
+ icorner_iglob(ncorner) = iglob
+ which_PML_poin(ibound,iglob) = .true.
+ end do
+ end do
+ nspec_PML=nspec_PML+1
+ end if
+ end do
+
+ end do !end nelem_thickness loop
+
+ write(IOUT,*) "number of PML spectral elements on side ", ibound,":", nspec_PML
+
+ enddo ! end loop on the 4 boundaries
+
+ !construction of table to use less memory for absorbing coefficients
+ ! allocate(spec_to_PML(nspec))
+ spec_to_PML=0
+ nspec_PML=0
+ do ispec=1,nspec
+ if (is_PML(ispec)) then
+ nspec_PML=nspec_PML+1
+ spec_to_PML(ispec)=nspec_PML
+ end if
+ enddo
+
+ ! allocate(ibool_PML(NGLLX,NGLLZ,nspec))
+ if (nspec_PML > 0) then
+ allocate(iPML_to_iglob(nspec_PML*NGLLX*NGLLZ))
+ else
+ allocate(iPML_to_iglob(1))
+ end if
+
+ iPML_to_iglob(:)=0
+ ibool_PML=0
+
+ npoin_PML=0
+ do ispec=1,nspec
+ if (is_PML(ispec)) then
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob=ibool(i,j,ispec)
+ k=1
+ do while (iglob/=iPML_to_iglob(k).and.iPML_to_iglob(k)/=0)
+ k=k+1
+ end do
+ ibool_PML(i,j,ispec)=k
+ if (k>npoin_PML) then
+ npoin_PML=npoin_PML+1
+ iPML_to_iglob(k)=iglob
+ end if
+ end do
+ end do
+ end if
+ end do
+
+ deallocate(iPML_to_iglob)
+ ! deallocate(icorner_iglob)
+ ! deallocate(which_PML_poin)
+
+ write(IOUT,*) "number of PML spectral elements :", nspec_PML
+ write(IOUT,*) "number of PML spectral points :", npoin_PML
+
+ end subroutine pml_init
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+!!!!!!!!!! subroutine define_PML_coefficients(npoin,nspec,is_PML,a_x,a_z,b_x,b_z,ibool,coord,&
+ subroutine define_PML_coefficients(npoin,nspec,is_PML,ibool,coord,&
+ which_PML_elem,kmato,density,poroelastcoef,numat,f0_temp,npoin_PML,&
+ ibool_PML,myrank,&
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store)
+
+ implicit none
+
+ include "constants.h"
+
+#ifdef USE_MPI
+ include "mpif.h"
+#endif
+
+ integer nspec, npoin, i, j,numat, ispec,iglob,npoin_PML,iPML
+ double precision :: f0_temp !xiezhinan
+
+ logical, dimension(nspec) :: is_PML
+ logical, dimension(4,nspec) :: which_PML_elem
+!!!!!!!!!!!! double precision, dimension(npoin_PML) :: a_x,a_z,b_x,b_z
+!!!!!!!!!!!!! DK DK added this
+ real(kind=CUSTOM_REAL), dimension(npoin_PML) :: &
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store !xiezhinan
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: coord
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool, ibool_PML
+ double precision, dimension(2,numat) :: density
+ double precision, dimension(4,3,numat) :: poroelastcoef
+ integer, dimension(nspec) :: kmato
+
+ double precision :: d_x, d_z, K_x, K_z, alpha_x, alpha_z
+! define an alias for y and z variable names (which are the same)
+ double precision :: d_y, alpha_y, K_y!!!!!!!!!!,aa,bb
+ double precision :: d0_z_bottom, d0_x_right, d0_z_top, d0_x_left
+ double precision :: abscissa_in_PML, abscissa_normalized
+
+ double precision :: thickness_PML_z_min_bottom,thickness_PML_z_max_bottom,&
+ thickness_PML_x_min_right,thickness_PML_x_max_right,&
+ thickness_PML_z_min_top,thickness_PML_z_max_top,&
+ thickness_PML_x_min_left,thickness_PML_x_max_left,&
+ thickness_PML_z_bottom,thickness_PML_x_right,&
+ thickness_PML_z_top,thickness_PML_x_left
+
+ double precision :: xmin, xmax, zmin, zmax, vpmax, xval, zval
+ double precision :: xoriginleft, xoriginright, zorigintop, zoriginbottom
+
+#ifdef USE_MPI
+! for MPI and partitioning
+ integer :: ier
+ integer :: nproc
+ integer :: iproc
+ character(len=256) :: prname
+
+ double precision :: thickness_PML_z_min_bottom_glob,thickness_PML_z_max_bottom_glob,&
+ thickness_PML_x_min_right_glob,thickness_PML_x_max_right_glob,&
+ thickness_PML_z_min_top_glob,thickness_PML_z_max_top_glob,&
+ thickness_PML_x_min_left_glob,thickness_PML_x_max_left_glob,&
+ thickness_PML_z_bottom_glob,thickness_PML_x_right_glob,&
+ thickness_PML_z_top_glob,thickness_PML_x_left_glob
+
+ double precision :: xmin_glob, xmax_glob, zmin_glob, zmax_glob, vpmax_glob, d0
+#endif
+
+ integer :: myrank
+
+ !PML fixed parameters
+
+! power to compute d0 profile
+ double precision, parameter :: NPOWER = 2.d0
+
+ double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
+
+ double precision :: ALPHA_MAX_PML
+
+ double precision :: Rcoef
+
+ double precision :: factorx, factorz
+
+! material properties of the elastic medium
+ double precision :: lambdalplus2mul_relaxed
+
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
+ Rcoef = 0.001d0
+
+ ALPHA_MAX_PML = 2.d0*PI*(f0_temp/2.d0) ! from Festa and Vilotte
+
+! check that NPOWER is okay
+ if(NPOWER < 1) stop 'NPOWER must be greater than 1'
+
+ ! get minimum and maximum values of mesh coordinates
+ xmin = minval(coord(1,:))
+ zmin = minval(coord(2,:))
+ xmax = maxval(coord(1,:))
+ zmax = maxval(coord(2,:))
+
+#ifdef USE_MPI
+ call MPI_ALLREDUCE (xmin, xmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (zmin, zmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (xmax, xmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (zmax, zmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ xmin = xmin_glob
+ zmin = zmin_glob
+ xmax = xmax_glob
+ zmax = zmax_glob
+#endif
+
+ thickness_PML_z_min_bottom=1.d30
+ thickness_PML_z_max_bottom=-1.d30
+
+ thickness_PML_x_min_right=1.d30
+ thickness_PML_x_max_right=-1.d30
+
+ thickness_PML_z_min_top=1.d30
+ thickness_PML_z_max_top=-1.d30
+
+ thickness_PML_x_min_left=1.d30
+ thickness_PML_x_max_left=-1.d30
+
+ do ispec=1,nspec
+ if (is_PML(ispec)) then
+ do j=1,NGLLZ
+ do i=1,NGLLX
+
+!!!bottom_case
+ if (which_PML_elem(IBOTTOM,ispec)) then
+ thickness_PML_z_max_bottom=max(coord(2,ibool(i,j,ispec)),thickness_PML_z_max_bottom)
+ thickness_PML_z_min_bottom=min(coord(2,ibool(i,j,ispec)),thickness_PML_z_min_bottom)
+ endif
+
+!!!right case
+ if (which_PML_elem(IRIGHT,ispec)) then
+ thickness_PML_x_max_right=max(coord(1,ibool(i,j,ispec)),thickness_PML_x_max_right)
+ thickness_PML_x_min_right=min(coord(1,ibool(i,j,ispec)),thickness_PML_x_min_right)
+ endif
+
+!!!top case
+ if (which_PML_elem(ITOP,ispec)) then
+ thickness_PML_z_max_top=max(coord(2,ibool(i,j,ispec)),thickness_PML_z_max_top)
+ thickness_PML_z_min_top=min(coord(2,ibool(i,j,ispec)),thickness_PML_z_min_top)
+ endif
+
+!!!left case
+ if (which_PML_elem(ILEFT,ispec)) then
+ thickness_PML_x_max_left=max(coord(1,ibool(i,j,ispec)),thickness_PML_x_max_left)
+ thickness_PML_x_min_left=min(coord(1,ibool(i,j,ispec)),thickness_PML_x_min_left)
+ endif
+
+ enddo
+ enddo
+ endif
+ enddo
+
+#ifdef USE_MPI
+
+!!!bottom
+ call MPI_ALLREDUCE (thickness_PML_z_max_bottom, thickness_PML_z_max_bottom_glob, &
+ 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (thickness_PML_z_min_bottom, thickness_PML_z_min_bottom_glob, &
+ 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ thickness_PML_z_max_bottom=thickness_PML_z_max_bottom_glob
+ thickness_PML_z_min_bottom=thickness_PML_z_min_bottom_glob
+
+!!!right
+ call MPI_ALLREDUCE (thickness_PML_x_max_right, thickness_PML_x_max_right_glob, &
+ 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (thickness_PML_x_min_right, thickness_PML_x_min_right_glob, &
+ 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ thickness_PML_x_max_right=thickness_PML_x_max_right_glob
+ thickness_PML_x_min_right=thickness_PML_x_min_right_glob
+
+!!!top
+ call MPI_ALLREDUCE (thickness_PML_z_max_top, thickness_PML_z_max_top_glob, &
+ 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (thickness_PML_z_min_top, thickness_PML_z_min_top_glob, &
+ 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ thickness_PML_z_max_top=thickness_PML_z_max_top_glob
+ thickness_PML_z_min_top=thickness_PML_z_min_top_glob
+
+!!!left
+ call MPI_ALLREDUCE (thickness_PML_x_max_left, thickness_PML_x_max_left_glob, &
+ 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (thickness_PML_x_min_left, thickness_PML_x_min_left_glob, &
+ 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ thickness_PML_x_max_left=thickness_PML_x_max_left_glob
+ thickness_PML_x_min_left=thickness_PML_x_min_left_glob
+#endif
+
+ thickness_PML_x_left= thickness_PML_x_max_left - thickness_PML_x_min_left
+ thickness_PML_x_right= thickness_PML_x_max_right - thickness_PML_x_min_right
+ thickness_PML_z_bottom= thickness_PML_z_max_bottom - thickness_PML_z_min_bottom
+ thickness_PML_z_top= thickness_PML_z_max_top - thickness_PML_z_min_top
+
+ vpmax = 0
+ do ispec = 1,nspec
+ if (is_PML(ispec)) then
+ ! get relaxed elastic parameters of current spectral element
+ lambdalplus2mul_relaxed = poroelastcoef(3,1,kmato(ispec))
+ vpmax=max(vpmax,sqrt(lambdalplus2mul_relaxed/density(1,kmato(ispec))))
+ endif
+ enddo
+
+#ifdef USE_MPI
+ call MPI_ALLREDUCE (vpmax, vpmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ vpmax=vpmax_glob
+#endif
+
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
+ d0_x_left = - (NPOWER + 1) * vpmax * log(Rcoef) / (2.d0 * thickness_PML_x_left)
+ d0_x_right = - (NPOWER + 1) * vpmax * log(Rcoef) / (2.d0 * thickness_PML_x_right)
+ d0_z_bottom = - (NPOWER + 1) * vpmax * log(Rcoef) / (2.d0 * thickness_PML_z_bottom)
+ d0_z_top = - (NPOWER + 1) * vpmax * log(Rcoef) / (2.d0 * thickness_PML_z_top)
+
+ if (myrank == 0) then
+ write(IOUT,*)
+ write(IOUT,*) 'PML properties -------'
+ write(IOUT,*) ' Vpmax=', vpmax
+ write(IOUT,*) ' log(Rcoef)=',log(Rcoef)
+ write(IOUT,*) ' thickness_PML_z_bottom =',thickness_PML_z_bottom
+ write(IOUT,*) ' thickness_PML_x_right =',thickness_PML_x_right
+ write(IOUT,*) ' thickness_PML_z_top =',thickness_PML_z_top
+ write(IOUT,*) ' thickness_PML_x_left =',thickness_PML_x_left
+ write(IOUT,*) ' d0_bottom =', d0_z_bottom
+ write(IOUT,*) ' d0_right =', d0_x_right
+ write(IOUT,*) ' d0_top =', d0_z_top
+ write(IOUT,*) ' d0_left =', d0_x_left
+ endif
+
+ d_x = ZERO
+ d_z = ZERO
+ K_x = ZERO
+ K_z = ZERO
+ alpha_x = ZERO
+ alpha_z = ZERO
+
+! origin of the PML layer (position of right edge minus thickness, in meters)
+ xoriginleft = thickness_PML_x_left+xmin
+ xoriginright = xmax - thickness_PML_x_right
+ zoriginbottom = thickness_PML_z_bottom + zmin
+ zorigintop = zmax-thickness_PML_z_top
+
+ do ispec = 1,nspec
+
+ if (is_PML(ispec)) then
+
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob=ibool(i,j,ispec)
+!! DK DK il y a un autre bug (un troisieme), c'est que certaines valeurs de ibool_PML() sont egales a zero
+!! DK DK ce qui n'est pas normal car ca doit etre seulement des numeros de points.
+!! DK DK cela n'est probablement pas grave car ca doit correspondre aux elements qui ne sont pas dans les PMLs, donc ignores
+ iPML=ibool_PML(i,j,ispec)
+!! DK DK added this
+ if(iPML < 1) stop 'iPML < 1 in a PML element'
+!! DK DK added this
+ ! abscissa of current grid point along the damping profile
+ xval = coord(1,iglob)
+ zval = coord(2,iglob)
+
+!!!! ---------- bottom edge
+ if (which_PML_elem(IBOTTOM,ispec)) then
+ ! define damping profile at the grid points
+ abscissa_in_PML = zoriginbottom - zval
+ if(abscissa_in_PML > 0.d0) then
+ abscissa_normalized = abscissa_in_PML / thickness_PML_z_bottom
+ d_z = d0_z_bottom / 0.6d0 * abscissa_normalized**NPOWER
+! alpha_z = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+! write(IOUT,*)d0_z_bottom,"d0_z_bottom"
+
+ alpha_z = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
+ + ALPHA_MAX_PML / 2.d0
+
+ K_z = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
+ else
+ d_z = 0.d0
+ alpha_z = 0.d0
+ K_z = 1.d0
+ endif
+ factorz = 1.d0
+ if (which_PML_elem(IRIGHT,ispec) .or. which_PML_elem(ILEFT,ispec)) then
+ factorx = 1.d0
+ else
+ factorx = 0.d0
+ endif
+ endif
+
+!!!! ---------- top edge
+ if (which_PML_elem(ITOP,ispec)) then
+ ! define damping profile at the grid points
+ abscissa_in_PML = zval - zorigintop
+ if(abscissa_in_PML > 0.d0) then
+ abscissa_normalized = abscissa_in_PML / thickness_PML_z_top
+ d_z = d0_z_top / 0.6d0 * abscissa_normalized**NPOWER
+! alpha_z = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+! write(IOUT,*)d0_z_top,"d0_z_top"
+
+ alpha_z = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
+ + ALPHA_MAX_PML / 2.d0
+
+ K_z = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
+ else
+ d_z = 0.d0
+ alpha_z = 0.d0
+ K_z = 1.d0
+ endif
+ factorz = 1.d0
+ if (which_PML_elem(IRIGHT,ispec) .or. which_PML_elem(ILEFT,ispec)) then
+ factorx = 1.d0
+ else
+ factorx = 0.d0
+ endif
+ endif
+
+!!!! ---------- right edge
+ if (which_PML_elem(IRIGHT,ispec)) then
+ ! define damping profile at the grid points
+ abscissa_in_PML = xval - xoriginright
+ if(abscissa_in_PML > 0.d0) then
+ abscissa_normalized = abscissa_in_PML / thickness_PML_x_right
+ d_x = d0_x_right / 0.6d0 * abscissa_normalized**NPOWER
+! write(IOUT,*)d0_x_right,"d0_x_right"
+! alpha_x = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+
+ alpha_x = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
+ + ALPHA_MAX_PML / 2.d0
+
+ K_x = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
+ else
+ d_x = 0.d0
+ alpha_x = 0.d0
+ K_x = 1.d0
+ endif
+ factorx = 1.d0
+ if (which_PML_elem(IBOTTOM,ispec) .or. which_PML_elem(ITOP,ispec)) then
+ factorz = 1.d0
+ else
+ factorz = 0.d0
+ endif
+ endif
+
+!!!! ---------- left edge
+ if (which_PML_elem(ILEFT,ispec)) then
+ ! define damping profile at the grid points
+ abscissa_in_PML = xoriginleft - xval
+ if(abscissa_in_PML > 0.d0) then
+ abscissa_normalized = abscissa_in_PML / thickness_PML_x_left
+ d_x = d0_x_left / 0.6d0 * abscissa_normalized**NPOWER
+! write(IOUT,*)d0_x_left,"d0_x_left"
+! alpha_x = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+
+ alpha_x = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
+ + ALPHA_MAX_PML / 2.d0
+
+ K_x = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
+ else
+ d_x = 0.d0
+ alpha_x = 0.d0
+ K_x = 1.d0
+ endif
+ factorx = 1.d0
+ if (which_PML_elem(IBOTTOM,ispec) .or. which_PML_elem(ITOP,ispec)) then
+ factorz = 1.d0
+ else
+ factorz = 0.d0
+ endif
+ endif
+
+!! DK DK define an alias for y and z variable names (which are the same)
+ K_y = K_z
+ d_y = d_z
+ alpha_y = alpha_z
+
+
+
+ K_x_store(iPML) = K_x
+ K_z_store(iPML) = K_y
+
+ d_x_store(iPML) = d_x
+ d_z_store(iPML) = d_y
+
+ alpha_x_store(iPML) = alpha_x
+ alpha_z_store(iPML) = alpha_y
+
+ write(IOUT,*)K_x_store(iPML),K_z_store(iPML),'K_store'
+ write(IOUT,*)d_x_store(iPML),d_z_store(iPML),'d_store'
+ write(IOUT,*)alpha_x_store(iPML),alpha_z_store(iPML),'alpha_store'
+
+ enddo
+ enddo
+ endif
+ enddo
+
+!!!!!!! write(IOUT,*) 'max bx',maxval(b_x),'min bx', minval(b_x)
+!!!!!!! write(IOUT,*) 'max ax',maxval(a_x),'min ax', minval(a_x)
+!!!!!!! write(IOUT,*) 'max by',maxval(b_z),'min by', minval(b_z)
+!!!!!!! write(IOUT,*) 'max ay',maxval(a_z),'min ay', minval(a_z)
+
+ end subroutine define_PML_coefficients
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-06-09 00:55:48 UTC (rev 20341)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-06-10 21:11:46 UTC (rev 20342)
@@ -482,8 +482,8 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_acoustic_adj_coupling
! inverse mass matrices
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic_one !zhinan
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic_three !zhinan
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic_one
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic_three
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_acoustic
real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
@@ -956,6 +956,30 @@
integer :: count_nspec_acoustic,count_nspec_acoustic_total,nspec_total,nglob_total,nb_acoustic_DOFs,nb_elastic_DOFs
double precision :: ratio_1DOF,ratio_2DOFs
+!PML parameters
+ logical, dimension(:), allocatable :: is_PML
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
+
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ rmemory_dux_dx,rmemory_duz_dx,rmemory_dux_dz,rmemory_duz_dz, &
+ rmemory_dux_dx_corner,rmemory_duz_dx_corner,rmemory_dux_dz_corner,rmemory_duz_dz_corner
+
+ integer, dimension(:), allocatable :: spec_to_PML,icorner_iglob
+ integer, dimension(:,:,:), allocatable :: ibool_PML
+ logical, dimension(:,:), allocatable :: which_PML_elem, which_PML_poin
+
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic_corner
+
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rmemory_potential_acoustic,&
+ rmemory_potential_ac_corner,&
+ rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
+ rmemory_acoustic_dux_dx_corner,rmemory_acoustic_dux_dz_corner
+
+ logical :: anyabs_glob
+ integer :: nspec_PML, npoin_PML
+
!***********************************************************************
!
! i n i t i a l i z a t i o n p h a s e
@@ -2742,6 +2766,165 @@
allocate(rhorho_ac_hessian_final1(1,1,1))
endif
+ ! PML absorbing conditionds
+ anyabs_glob=anyabs
+#ifdef USE_MPI
+ call MPI_ALLREDUCE(anyabs, anyabs_glob, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ier)
+#endif
+
+ if( PML_BOUNDARY_CONDITIONS .and. anyabs_glob ) then
+
+ !PML code
+ allocate(is_PML(nspec))
+ allocate(icorner_iglob(nglob))
+ allocate(which_PML_elem(4,nspec))
+ allocate(which_PML_poin(4,nglob))
+ allocate(spec_to_PML(nspec))
+ allocate(ibool_PML(NGLLX,NGLLZ,nspec))
+
+ is_PML(:) = .false.
+ which_PML_elem(:,:) = .false.
+ which_PML_poin(:,:) = .false.
+
+ call pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
+ nspec_PML,is_PML,which_PML_elem,which_PML_poin,spec_to_PML,ibool_PML, &
+ npoin_PML,icorner_iglob)
+
+ deallocate(icorner_iglob)
+ deallocate(which_PML_poin)
+
+ if (npoin_PML==0) npoin_PML=1
+!!!!!!!!!!!!! DK DK added this
+
+ if (npoin_PML > 0) then
+
+ allocate(K_x_store(npoin_PML))
+ allocate(K_z_store(npoin_PML))
+ allocate(d_x_store(npoin_PML))
+ allocate(d_z_store(npoin_PML))
+ allocate(alpha_x_store(npoin_PML))
+ allocate(alpha_z_store(npoin_PML))
+
+ !! DK DK initialize to zero
+ K_x_store(:) = 0
+ K_z_store(:) = 0
+ d_x_store(:) = 0
+ d_z_store(:) = 0
+ alpha_x_store(:) = 0
+ alpha_z_store(:) = 0
+
+ call define_PML_coefficients(nglob,nspec,is_PML,ibool,coord,&
+ which_PML_elem,kmato,density,poroelastcoef,numat,f0(1),npoin_PML,&
+ ibool_PML,myrank,&
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store)
+
+ endif
+
+ !elastic PML memory variables
+ if (any_elastic .and. nspec_PML>0) then
+
+ allocate(rmemory_displ_elastic(2,3,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_displ_elastic_corner(2,3,NGLLX,NGLLZ,nspec_PML))
+
+ allocate(rmemory_dux_dx(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_dux_dz(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_duz_dx(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_duz_dz(2,NGLLX,NGLLZ,nspec_PML))
+
+ allocate(rmemory_dux_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_dux_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_duz_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_duz_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
+
+ else
+
+ allocate(rmemory_displ_elastic(1,1,1,1,1))
+ allocate(rmemory_displ_elastic_corner(1,1,1,1,1))
+
+ allocate(rmemory_dux_dx(1,1,1,1))
+ allocate(rmemory_dux_dz(1,1,1,1))
+ allocate(rmemory_duz_dx(1,1,1,1))
+ allocate(rmemory_duz_dz(1,1,1,1))
+
+ allocate(rmemory_dux_dx_corner(1,1,1,1))
+ allocate(rmemory_dux_dz_corner(1,1,1,1))
+ allocate(rmemory_duz_dx_corner(1,1,1,1))
+ allocate(rmemory_duz_dz_corner(1,1,1,1))
+
+ end if
+
+ if (any_acoustic .and. nspec_PML>0) then
+
+ allocate(rmemory_potential_acoustic(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_potential_ac_corner(2,NGLLX,NGLLZ,nspec_PML))
+
+ allocate(rmemory_acoustic_dux_dx(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_acoustic_dux_dz(2,NGLLX,NGLLZ,nspec_PML))
+
+ allocate(rmemory_acoustic_dux_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_acoustic_dux_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
+
+ else
+
+ allocate(rmemory_potential_acoustic(1,1,1,1))
+ allocate(rmemory_potential_ac_corner(1,1,1,1))
+
+ allocate(rmemory_acoustic_dux_dx(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dz(1,1,1,1))
+
+ allocate(rmemory_acoustic_dux_dx_corner(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dz_corner(1,1,1,1))
+
+ end if
+
+ else
+
+
+!!!!!!!!!!!!! DK DK added this
+
+ allocate(rmemory_dux_dx(1,1,1,1))
+ allocate(rmemory_dux_dz(1,1,1,1))
+ allocate(rmemory_duz_dx(1,1,1,1))
+ allocate(rmemory_duz_dz(1,1,1,1))
+
+ allocate(rmemory_dux_dx_corner(1,1,1,1))
+ allocate(rmemory_dux_dz_corner(1,1,1,1))
+ allocate(rmemory_duz_dx_corner(1,1,1,1))
+ allocate(rmemory_duz_dz_corner(1,1,1,1))
+
+ allocate(rmemory_displ_elastic(1,1,1,1,1))
+ allocate(rmemory_displ_elastic_corner(1,1,1,1,1))
+
+ allocate(rmemory_acoustic_dux_dx(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dz(1,1,1,1))
+
+ allocate(rmemory_acoustic_dux_dx_corner(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dz_corner(1,1,1,1))
+
+ allocate(is_PML(1))
+ allocate(ibool_PML(1,1,1))
+ allocate(spec_to_PML(1))
+ allocate(which_PML_elem(1,1))
+
+ allocate(K_x_store(1))
+ allocate(K_z_store(1))
+ allocate(d_x_store(1))
+ allocate(d_z_store(1))
+ allocate(alpha_x_store(1))
+ allocate(alpha_z_store(1))
+
+ end if ! PML_BOUNDARY_CONDITIONS
+
+ rmemory_dux_dx(:,:,:,:) = ZERO
+ rmemory_dux_dz(:,:,:,:) = ZERO
+ rmemory_duz_dx(:,:,:,:) = ZERO
+ rmemory_duz_dz(:,:,:,:) = ZERO
+
+ rmemory_dux_dx_corner(:,:,:,:) = ZERO
+ rmemory_dux_dz_corner(:,:,:,:) = ZERO
+ rmemory_duz_dx_corner(:,:,:,:) = ZERO
+ rmemory_duz_dz_corner(:,:,:,:) = ZERO
+
endif ! ipass == 1
!
@@ -5165,7 +5348,12 @@
nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k, &
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
- stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring)
+ stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring, &
+ is_PML,nspec_PML,npoin_PML,ibool_PML,spec_to_PML,which_PML_elem, &
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
+ rmemory_displ_elastic,rmemory_displ_elastic_corner,&
+ rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz,&
+ rmemory_dux_dx_corner,rmemory_dux_dz_corner,rmemory_duz_dx_corner,rmemory_duz_dz_corner)
if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
!--- left absorbing boundary
More information about the CIG-COMMITS
mailing list