[cig-commits] r20481 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Sat Jul 7 04:56:24 PDT 2012
Author: xie.zhinan
Date: 2012-07-07 04:56:23 -0700 (Sat, 07 Jul 2012)
New Revision: 20481
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
add pml for fluid and fluid-solid simulation. add some flags to stop using PML in poroelastic case and in mpiversion of the code.
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-07-07 02:09:18 UTC (rev 20480)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-07-07 11:56:23 UTC (rev 20481)
@@ -54,7 +54,13 @@
SIMULATION_TYPE,SAVE_FORWARD,nspec_left,nspec_right,&
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
- b_absorb_acoustic_bottom,b_absorb_acoustic_top,IS_BACKWARD_FIELD)
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top,IS_BACKWARD_FIELD,&
+ 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_potential_acoustic,rmemory_potential_acoustic_corner,&
+ rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
+ rmemory_acoustic_dux_dx_corner,rmemory_acoustic_dux_dz_corner,deltat,&
+ PML_BOUNDARY_CONDITIONS)
! compute forces for the acoustic elements
@@ -121,6 +127,28 @@
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,NGLLX,NGLLZ,nspec_PML) :: rmemory_potential_acoustic,rmemory_potential_acoustic_corner
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLZ,nspec_PML) :: &
+ rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
+ rmemory_acoustic_dux_dx_corner,rmemory_acoustic_dux_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(NGLLX,NGLLZ,nspec_PML) :: potential_dot_dot_acoustic_PML, &
+ potential_dot_dot_acoustic_PML_corner
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) ::PML_dux_dxl,PML_dux_dzl,&
+ PML_dux_dxl_new,PML_dux_dzl_new
+ real(kind=CUSTOM_REAL) :: coef0_x, coef1_x, coef2_x, coef0_z, coef1_z, coef2_z,bb,aa,&
+ deltat
+
+ logical :: PML_BOUNDARY_CONDITIONS
+
ifirstelem = 1
ilastelem = nspec
@@ -157,6 +185,278 @@
! derivatives of potential
dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+
+
+ ! derivative along x and along z
+ if(PML_BOUNDARY_CONDITIONS .and. is_PML(ispec))then
+
+ ispec_PML=spec_to_PML(ispec)
+
+ PML_dux_dxl(i,j,ispec_PML) = dux_dxl
+ PML_dux_dzl(i,j,ispec_PML)=dux_dzl
+
+ dux_dxi = ZERO
+ dux_dgamma = 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 = dux_dxi &
+ +(potential_acoustic(ibool(k,j,ispec))+deltat*potential_dot_acoustic(ibool(k,j,ispec)))*hprime_xx(i,k)
+ dux_dgamma = dux_dgamma &
+ +(potential_acoustic(ibool(i,k,ispec))+deltat*potential_dot_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)
+
+ ! derivatives of potential
+ dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
+ dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+
+ PML_dux_dxl_new(i,j,ispec_PML) = dux_dxl
+ PML_dux_dzl_new(i,j,ispec_PML) = dux_dzl
+ endif
+
+
+ if(PML_BOUNDARY_CONDITIONS .and. is_PML(ispec))then
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ if (which_PML_elem(ILEFT,ispec) .or. which_PML_elem(IRIGHT,ispec)) then
+
+ bb = 0.d0 - (d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML))
+ if(abs(alpha_x_store(iPML)) > 1.d-3) then
+ aa = - d_x_store(iPML) / &
+ (k_x_store(iPML) * d_x_store(iPML) + k_x_store(iPML)**2 * alpha_x_store(iPML))
+ else
+ aa = - 1.d0 / k_x_store(iPML)
+ endif
+
+ coef0_x = exp(bb * deltat)
+ coef1_x = (1.d0 - exp(bb * deltat / 2.d0)) * aa
+ coef2_x = (1.d0 - exp(bb* deltat / 2.d0)) * exp(bb * deltat / 2.d0) * aa
+
+
+ rmemory_acoustic_dux_dx(1,i,j,ispec_PML) = 0.d0
+ rmemory_acoustic_dux_dx(2,i,j,ispec_PML) = coef0_x*rmemory_acoustic_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
+
+ if(abs(alpha_x_store(iPML)) > 1.d-3) then
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * d_x_store(iPML)/bb
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * &
+ d_x_store(iPML)/bb
+ else
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = deltat / 2.0d0 * d_x_store(iPML)
+ coef2_x = deltat / 2.0d0 * d_x_store(iPML)
+ endif
+
+ rmemory_acoustic_dux_dz(1,i,j,ispec_PML) = coef0_x * rmemory_acoustic_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_acoustic_dux_dz(2,i,j,ispec_PML) = 0.d0
+
+ dux_dxl = PML_dux_dxl(i,j,ispec_PML) + rmemory_acoustic_dux_dx(1,i,j,ispec_PML) &
+ + rmemory_acoustic_dux_dx(2,i,j,ispec_PML)
+ dux_dzl = PML_dux_dzl(i,j,ispec_PML) + rmemory_acoustic_dux_dz(1,i,j,ispec_PML) &
+ + rmemory_acoustic_dux_dz(2,i,j,ispec_PML)
+
+ 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)
+
+ if(abs(alpha_z_store(iPML)) .lt. 1.d-3)then
+ bb = alpha_z_store(iPML)
+ if(abs(k_x_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))+d_x_store(iPML)) .lt. 1.d-3)then
+ aa=d_z_store(iPML)/k_x_store(iPML)
+ else
+ aa = 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))
+ endif
+ coef0_x = exp(- bb * deltat)
+ coef1_x = deltat / 2.0d0 * aa
+ coef2_x = deltat / 2.0d0 * aa
+ else
+ if(abs(k_x_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))+d_x_store(iPML)) .lt. 1.d-3)then
+ aa=d_z_store(iPML)/k_x_store(iPML)
+ else
+ aa = 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))
+ endif
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * aa/bb
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * aa /bb
+ endif
+
+
+ rmemory_acoustic_dux_dx_corner(1,i,j,ispec_PML) = coef0_x*rmemory_acoustic_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
+
+
+ if(abs(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)*(alpha_x_store(iPML)-alpha_z_store(iPML))+d_x_store(iPML))/&
+ k_x_store(iPML)**2) .lt. 1.d-3 .and. &
+ abs(k_x_store(iPML)*alpha_x_store(iPML)+d_x_store(iPML)) .lt. 1.d-3)then
+ coef0_z = 1.d0
+ coef1_z = 0.d0
+ coef2_z = 0.d0
+
+ rmemory_acoustic_dux_dx_corner(2,i,j,ispec_PML) = coef0_z*rmemory_acoustic_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
+
+
+ elseif(abs(k_x_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))+d_x_store(iPML)) .lt. 1.d-3 .and. &
+ abs(k_x_store(iPML)*alpha_x_store(iPML)+d_x_store(iPML)) .lt. 1.d-3)then
+ bb=(k_x_store(iPML)*alpha_x_store(iPML)+d_x_store(iPML))/k_x_store(iPML)
+ aa=d_z_store(iPML)/k_x_store(iPML)
+ coef0_z = 1.d0
+ coef1_z = 0.5d0 *deltat
+ coef2_z = 0.5d0 *deltat
+
+ rmemory_acoustic_dux_dx_corner(2,i,j,ispec_PML) = coef0_z*rmemory_acoustic_dux_dx_corner(2,i,j,ispec_PML) &
+ + PML_dux_dxl_new(i,j,ispec_PML) * coef1_z*aa + PML_dux_dxl(i,j,ispec_PML) * coef2_z*aa
+
+
+ else
+ bb=(k_x_store(iPML)*alpha_x_store(iPML)+d_x_store(iPML))/k_x_store(iPML)
+ aa=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)*(alpha_x_store(iPML)-alpha_z_store(iPML))+d_x_store(iPML))/&
+ k_x_store(iPML)**2
+ coef0_z = exp(- bb * deltat)
+ coef1_z = (1 - exp(- bb * deltat / 2)) * aa/bb
+ coef2_z = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * aa /bb
+
+ rmemory_acoustic_dux_dx_corner(2,i,j,ispec_PML) = coef0_z*rmemory_acoustic_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
+
+
+ endif
+
+ if(abs(alpha_x_store(iPML)) .lt. 1.d-3)then
+ bb = alpha_x_store(iPML)
+ if(abs(k_z_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))-d_z_store(iPML)) .lt. 1.d-3)then
+ aa=d_x_store(iPML)/k_z_store(iPML)
+ else
+ aa = d_x_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))/&
+ (k_z_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))-d_z_store(iPML))
+ endif
+ coef0_x = exp(- bb * deltat)
+ coef1_x = deltat / 2.0d0 * aa
+ coef2_x = deltat / 2.0d0 * aa
+ else
+ if(abs(k_z_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))-d_z_store(iPML)) .lt. 1.d-3)then
+ aa=d_x_store(iPML)/k_z_store(iPML)
+ else
+ aa = d_x_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))/&
+ (k_z_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))-d_z_store(iPML))
+ endif
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * aa/bb
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * aa /bb
+ endif
+
+ rmemory_acoustic_dux_dz_corner(1,i,j,ispec_PML) = coef0_x * rmemory_acoustic_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
+
+
+ if(abs(d_z_store(iPML)*&
+ (k_x_store(iPML)*k_z_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)*(alpha_z_store(iPML)-alpha_x_store(iPML))+d_z_store(iPML))/&
+ k_z_store(iPML)**2) .lt. 1.d-3 .and. &
+ abs(k_z_store(iPML)*alpha_z_store(iPML)+d_z_store(iPML)) .lt. 1.d-3)then
+ coef0_z = 1.d0
+ coef1_z = 0.d0
+ coef2_z = 0.d0
+
+ rmemory_acoustic_dux_dz_corner(2,i,j,ispec_PML) = coef0_z * rmemory_acoustic_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
+
+ elseif(abs(k_z_store(iPML)*(alpha_z_store(iPML)-alpha_x_store(iPML))+d_z_store(iPML)) .lt. 1.d-3 .and. &
+ abs(k_z_store(iPML)*alpha_z_store(iPML)+d_z_store(iPML)) .lt. 1.d-3)then
+ bb=(k_z_store(iPML)*alpha_z_store(iPML)+d_z_store(iPML))/k_z_store(iPML)
+ aa=d_x_store(iPML)/k_z_store(iPML)
+ coef0_z = 1.d0
+ coef1_z = 0.5d0 *deltat
+ coef2_z = 0.5d0 *deltat
+
+ rmemory_acoustic_dux_dz_corner(2,i,j,ispec_PML) = coef0_z * rmemory_acoustic_dux_dz_corner(2,i,j,ispec_PML) &
+ + PML_dux_dzl_new(i,j,ispec_PML) *coef1_z*aa + PML_dux_dzl(i,j,ispec_PML) * coef2_z*aa
+
+ else
+ bb=(k_z_store(iPML)*alpha_z_store(iPML)+d_z_store(iPML))/k_z_store(iPML)
+ aa=d_z_store(iPML)*&
+ (k_x_store(iPML)*k_z_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)*(alpha_z_store(iPML)-alpha_x_store(iPML))+d_z_store(iPML))/&
+ k_z_store(iPML)**2
+ coef0_z = exp(- bb * deltat)
+ coef1_z = (1 - exp(- bb * deltat / 2)) * aa/bb
+ coef2_z = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * aa /bb
+
+ rmemory_acoustic_dux_dz_corner(2,i,j,ispec_PML) = coef0_z * rmemory_acoustic_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
+
+ endif
+
+ dux_dxl = PML_dux_dxl(i,j,ispec_PML) + rmemory_acoustic_dux_dx_corner(1,i,j,ispec_PML) +&
+ rmemory_acoustic_dux_dx_corner(2,i,j,ispec_PML)
+ dux_dzl = PML_dux_dzl(i,j,ispec_PML) + rmemory_acoustic_dux_dz_corner(1,i,j,ispec_PML) + &
+ rmemory_acoustic_dux_dz_corner(2,i,j,ispec_PML)
+
+ endif
+
+ else
+
+ if(abs(alpha_z_store(iPML)) > 1.d-3) then
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * d_z_store(iPML)/bb
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * &
+ d_z_store(iPML) / bb
+ else
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = deltat / 2.0d0 * d_z_store(iPML)
+ coef2_x = deltat / 2.0d0 * d_z_store(iPML)
+ endif
+
+
+ rmemory_acoustic_dux_dx(1,i,j,ispec_PML) = coef0_x*rmemory_acoustic_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_acoustic_dux_dx(2,i,j,ispec_PML) = 0.d0
+
+ bb = 0.d0 - (d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML))
+ if(abs(alpha_z_store(iPML)) > 1.d-3) then
+ aa = - d_z_store(iPML) * k_x_store(iPML) / &
+ (k_z_store(iPML) * d_z_store(iPML) + k_z_store(iPML)**2 * alpha_z_store(iPML))
+ else
+ aa = - k_x_store(iPML) / k_z_store(iPML)
+ endif
+ coef0_x = exp(bb * deltat)
+ coef1_x = (1.d0 - exp(bb * deltat / 2.d0)) * aa
+ coef2_x = (1.d0 - exp(bb* deltat / 2.d0)) * exp(bb * deltat / 2.d0) * aa
+
+ rmemory_acoustic_dux_dz(1,i,j,ispec_PML) = 0.d0
+ rmemory_acoustic_dux_dz(2,i,j,ispec_PML) = coef0_x * rmemory_acoustic_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
+
+ dux_dxl = dux_dxl + rmemory_acoustic_dux_dx(1,i,j,ispec_PML) + rmemory_acoustic_dux_dx(2,i,j,ispec_PML)
+ dux_dzl = dux_dzl + rmemory_acoustic_dux_dz(1,i,j,ispec_PML) + rmemory_acoustic_dux_dz(2,i,j,ispec_PML)
+
+ endif
+ endif
+
jacobianl = jacobian(i,j,ispec)
! if external density model
@@ -169,6 +469,159 @@
enddo
enddo
+
+ ! 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))
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec))
+ kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
+ rhol = density(1,kmato(ispec))
+ endif
+
+ if(is_PML(ispec))then
+ iPML=ibool_PML(i,j,ispec)
+ iglob=ibool(i,j,ispec)
+
+
+ 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)
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * (k_z_store(iPML)*d_x_store(iPML)*(bb))
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * &
+ (k_z_store(iPML)*d_x_store(iPML)*(bb))
+
+ rmemory_potential_acoustic(1,i,j,ispec_PML)=coef0_x * rmemory_potential_acoustic(1,i,j,ispec_PML) &
+ + (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob)) * coef1_x &
+ + potential_acoustic(iglob) * coef2_x
+
+ rmemory_potential_acoustic(2,i,j,ispec_PML)=0.0
+
+ 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)
+
+ if(abs(alpha_x_store(iPML)-alpha_z_store(iPML)).lt. 1.d-3)then
+ bb = alpha_x_store(iPML)
+ aa = (k_z_store(iPML)*d_x_store(iPML)*(bb))
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * aa
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * aa
+
+ rmemory_potential_acoustic_corner(1,i,j,ispec_PML)=&
+ coef0_x * rmemory_potential_acoustic_corner(1,i,j,ispec_PML) &
+ + (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob)) * coef1_x &
+ + potential_acoustic(iglob) * coef2_x
+ else
+ bb = alpha_x_store(iPML)
+ aa = bb*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))
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * aa
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * aa
+
+ rmemory_potential_acoustic_corner(1,i,j,ispec_PML)=&
+ coef0_x * rmemory_potential_acoustic_corner(1,i,j,ispec_PML) &
+ + (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob)) * coef1_x &
+ + potential_acoustic(iglob) * coef2_x
+
+ endif
+
+ if(abs(alpha_x_store(iPML)-alpha_z_store(iPML)).lt. 1.d-3)then
+ bb = alpha_z_store(iPML)
+ aa = k_x_store(iPML)*d_z_store(iPML)*(bb)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * aa
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * aa
+
+ rmemory_potential_acoustic_corner(2,i,j,ispec_PML)=&
+ coef0_x * rmemory_potential_acoustic_corner(2,i,j,ispec_PML) &
+ + (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob)) * coef1_x &
+ + potential_acoustic(iglob) * coef2_x
+
+ else
+ bb = alpha_z_store(iPML)
+ aa = bb*d_z_store(iPML)*(k_x_store(iPML)*(alpha_x_store(iPML)-alpha_z_store(iPML))+d_x_store(iPML))/&
+ (alpha_x_store(iPML)-alpha_z_store(iPML))
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * aa
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * aa
+
+ rmemory_potential_acoustic_corner(2,i,j,ispec_PML)=&
+ coef0_x * rmemory_potential_acoustic_corner(2,i,j,ispec_PML) &
+ + (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob)) * coef1_x &
+ + potential_acoustic(iglob) * coef2_x
+ endif
+
+
+ endif
+
+ else
+
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ iglob=ibool(i,j,ispec)
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ coef1_x = (1 - exp(- bb * deltat / 2)) * (k_x_store(iPML)*d_z_store(iPML)*(bb))
+ coef2_x = (1 - exp(- bb * deltat / 2)) * exp(- bb * deltat / 2) * &
+ (k_x_store(iPML)*d_z_store(iPML)*(bb))
+
+
+ rmemory_potential_acoustic(1,i,j,ispec_PML)=0.d0
+ rmemory_potential_acoustic(2,i,j,ispec_PML)=coef0_x * rmemory_potential_acoustic(2,i,j,ispec_PML) &
+ + (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob)) * coef1_x &
+ + potential_acoustic(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)
+ potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
+ ((d_x_store(iPML)*k_z_store(iPML))*potential_dot_acoustic(iglob)+ &
+ rmemory_potential_acoustic(1,i,j,ispec_PML)+rmemory_potential_acoustic(2,i,j,ispec_PML)&
+ -k_z_store(iPML)*d_x_store(iPML)*alpha_x_store(iPML)*potential_acoustic(iglob))
+
+
+ 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)
+
+ potential_dot_dot_acoustic_PML_corner(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
+ ((d_x_store(iPML)*k_z_store(iPML)+d_z_store(iPML)*k_x_store(iPML))*potential_dot_acoustic(iglob)+ &
+ rmemory_potential_acoustic_corner(1,i,j,ispec_PML)+rmemory_potential_acoustic_corner(2,i,j,ispec_PML)&
+ -k_z_store(iPML)*d_x_store(iPML)*alpha_x_store(iPML)*potential_acoustic(iglob)-&
+ k_x_store(iPML)*d_z_store(iPML)*alpha_z_store(iPML)*potential_acoustic(iglob)+&
+ d_z_store(iPML)*d_x_store(iPML)*potential_acoustic(iglob))
+
+ endif
+
+ else
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ iglob=ibool(i,j,ispec)
+ potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
+ ((d_z_store(iPML)*k_x_store(iPML))*potential_dot_acoustic(iglob)+ &
+ rmemory_potential_acoustic(1,i,j,ispec_PML)+rmemory_potential_acoustic(2,i,j,ispec_PML)&
+ -k_x_store(iPML)*d_z_store(iPML)*alpha_z_store(iPML)*potential_acoustic(iglob))
+ endif
+ endif
+
+ enddo
+ enddo
+
+
+
!
! second double-loop over GLL to compute all the terms
!
@@ -184,6 +637,32 @@
(tempx1(k,j)*hprimewgll_xx(k,i) + tempx2(i,k)*hprimewgll_zz(k,j))
enddo
+ 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)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - potential_dot_dot_acoustic_PML(i,j,ispec_PML)
+
+ 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)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ + potential_dot_dot_acoustic_PML(i,j,ispec_PML)
+
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - potential_dot_dot_acoustic_PML_corner(i,j,ispec_PML)
+ endif
+
+ else
+ ispec_PML=spec_to_PML(ispec)
+ iPML=ibool_PML(i,j,ispec)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - potential_dot_dot_acoustic_PML(i,j,ispec_PML)
+
+ endif
+ endif
+
enddo ! second loop over the GLL points
enddo
@@ -191,10 +670,71 @@
enddo ! end of loop over all spectral elements
+ if(PML_BOUNDARY_CONDITIONS .and. anyabs) 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)
+ potential_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_dot_acoustic(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)
+ potential_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_dot_acoustic(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)
+ potential_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_dot_acoustic(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)
+ potential_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_dot_dot_acoustic(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
+
+
!
!--- absorbing boundaries
!
- if(anyabs) then
+ if(.not. PML_BOUNDARY_CONDITIONS .and. anyabs) then
do ispecabs=1,nelemabs
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-07-07 02:09:18 UTC (rev 20480)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-07-07 11:56:23 UTC (rev 20481)
@@ -970,11 +970,9 @@
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 &
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 &
) then
- A3 = 0.d0
+ A3 = alpha_x_store(iPML) ** 2 * d_x_store(iPML) * k_z_store(iPML)
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) ) &
@@ -990,7 +988,6 @@
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
@@ -1008,11 +1005,9 @@
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 &
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 &
) then
- A4 = 0.d0
+ A4 = alpha_z_store(iPML) ** 2 * d_z_store(iPML) * k_x_store(iPML)
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) ) &
@@ -1740,7 +1735,7 @@
!
!--- absorbing boundaries
!
- if( PML_BOUNDARY_CONDITIONS ) then
+ if( PML_BOUNDARY_CONDITIONS .and. anyabs) then
! we have to put Dirichlet on the boundary of the PML
do ispecabs = 1,nelemabs
@@ -1748,7 +1743,7 @@
ispec = numabs(ispecabs)
if (is_PML(ispec)) then
- ispec_PML=spec_to_PML(ispec)
+ ispec_PML=spec_to_PML(ispec)
!--- left absorbing boundary
if(codeabs(ILEFT,ispecabs)) then
i = 1
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-07-07 02:09:18 UTC (rev 20480)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-07-07 11:56:23 UTC (rev 20481)
@@ -60,7 +60,7 @@
,coord &
#endif
,K_x_store,K_z_store,npoin_PML,ibool_PML,is_PML,&
- d_x_store,d_z_store,alpha_x_store,alpha_z_store,PML_BOUNDARY_CONDITIONS)
+ d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS)
! builds the global mass matrix
@@ -124,7 +124,7 @@
!!!!!!!!!!!!! DK DK added this
integer :: npoin_PML,iPML
- 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(npoin_PML) :: K_x_store,K_z_store,d_x_store,d_z_store
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool_PML
logical, dimension(nspec) :: is_PML
logical :: PML_BOUNDARY_CONDITIONS, anyabs_local
@@ -232,8 +232,18 @@
! for acoustic medium
+! rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+! + wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / kappal
+
+ if (is_PML(ispec)) then
+ iPML=ibool_PML(i,j,ispec)
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(iPML) * K_z_store(iPML)&
+ + (d_x_store(iPML)*k_z_store(iPML)+d_z_store(iPML)*k_x_store(iPML)) * deltat / 2.d0)
+ else
rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / kappal
+ endif
endif
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-07-07 02:09:18 UTC (rev 20480)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-07-07 11:56:23 UTC (rev 20481)
@@ -222,12 +222,12 @@
#endif
integer nspec, npoin, i, j,numat, ispec,iglob,npoin_PML,iPML
- double precision :: f0_temp !xiezhinan
+ double precision :: f0_temp
logical, dimension(nspec) :: is_PML
logical, dimension(4,nspec) :: which_PML_elem
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
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: coord
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool, ibool_PML
@@ -237,7 +237,6 @@
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
double precision :: d0_z_bottom, d0_x_right, d0_z_top, d0_x_left
double precision :: abscissa_in_PML, abscissa_normalized
@@ -254,16 +253,16 @@
#ifdef USE_MPI
! for MPI and partitioning
integer :: ier
- integer :: nproc
- integer :: iproc
- character(len=256) :: prname
+! 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
+ 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
@@ -445,6 +444,7 @@
zoriginbottom = thickness_PML_z_bottom + zmin
zorigintop = zmax-thickness_PML_z_top
+
do ispec = 1,nspec
if (is_PML(ispec)) then
@@ -462,11 +462,9 @@
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
+ 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
@@ -489,12 +487,9 @@
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
+ 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
@@ -516,12 +511,9 @@
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
+ 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
@@ -543,12 +535,9 @@
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
+ 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
@@ -566,33 +555,17 @@
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
-
+ K_z_store(iPML) = K_z
d_x_store(iPML) = d_x
- d_z_store(iPML) = d_y
-
+ d_z_store(iPML) = d_z
alpha_x_store(iPML) = alpha_x
- alpha_z_store(iPML) = alpha_y
+ alpha_z_store(iPML) = alpha_z
-! 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-07-07 02:09:18 UTC (rev 20480)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-07-07 11:56:23 UTC (rev 20481)
@@ -973,7 +973,7 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic_corner
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rmemory_potential_acoustic,&
- rmemory_potential_ac_corner,&
+ rmemory_potential_acoustic_corner,&
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
rmemory_acoustic_dux_dx_corner,rmemory_acoustic_dux_dz_corner
@@ -1018,11 +1018,18 @@
if(SIMULATION_TYPE == 2 .and.(time_stepping_scheme == 2 .or. time_stepping_scheme == 3)) &
stop 'RK and LDDRK time scheme not supported for adjoint inversion'
if(nproc /= nproc_read_from_database) stop 'must always have nproc == nproc_read_from_database'
+
+#ifdef USE_MPI
+ if(PML_BOUNDARY_CONDITIONS)then
+ stop 'PML_BOUNDARY_CONDITIONS do not ready for mpi version of the code'
+ endif
+#endif
+
!! DK DK Dec 2011: add a small crack (discontinuity) in the medium manually
npgeo_ori = npgeo
if(ADD_A_SMALL_CRACK_IN_THE_MEDIUM) npgeo = npgeo + NB_POINTS_TO_ADD_TO_NPGEO
-#ifndef USE_MPI
+#ifdef USE_MPI
if(PERFORM_CUTHILL_MCKEE) then
NUMBER_OF_PASSES = 2
else
@@ -1278,6 +1285,10 @@
anisotropic,elastic,poroelastic,porosity,anisotropy,kmato,numat, &
nspec,nspec_allocate,p_sv,ATTENUATION_VISCOELASTIC_SOLID,count_nspec_acoustic)
+ if(PML_BOUNDARY_CONDITIONS .and. any_poroelastic) then
+ stop 'PML boundary conditions do not ready for poroelastic simulation'
+ endif
+
! allocate memory variables for attenuation
if(ipass == 1) then
allocate(e1(NGLLX,NGLLZ,nspec_allocate,N_SLS))
@@ -2148,8 +2159,9 @@
xi_receiver,gamma_receiver,station_name,network_name,x_final_receiver,z_final_receiver,iglob_source)
! compute source array for adjoint source
+ if(ipass == 1) nadj_rec_local = 0
if(SIMULATION_TYPE == 2) then ! adjoint calculation
- nadj_rec_local = 0
+
do irec = 1,nrec
if(myrank == which_proc_receiver(irec)) then
! check that the source proc number is okay
@@ -2868,7 +2880,7 @@
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_potential_acoustic_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))
@@ -2879,7 +2891,7 @@
else
allocate(rmemory_potential_acoustic(1,1,1,1))
- allocate(rmemory_potential_ac_corner(1,1,1,1))
+ allocate(rmemory_potential_acoustic_corner(1,1,1,1))
allocate(rmemory_acoustic_dux_dx(1,1,1,1))
allocate(rmemory_acoustic_dux_dz(1,1,1,1))
@@ -4869,7 +4881,13 @@
SIMULATION_TYPE,SAVE_FORWARD,nspec_left,nspec_right,&
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
- b_absorb_acoustic_bottom,b_absorb_acoustic_top,.false.)
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top,.false.,&
+ 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_potential_acoustic,rmemory_potential_acoustic_corner,&
+ rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
+ rmemory_acoustic_dux_dx_corner,rmemory_acoustic_dux_dz_corner,deltat,&
+ PML_BOUNDARY_CONDITIONS)
if( SIMULATION_TYPE == 2 ) then
call compute_forces_acoustic(nglob,nspec,nelemabs,numat,it,NSTEP, &
anyabs,assign_external_model,ibool,kmato,numabs, &
@@ -4883,7 +4901,13 @@
SIMULATION_TYPE,SAVE_FORWARD,nspec_left,nspec_right,&
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
- b_absorb_acoustic_bottom,b_absorb_acoustic_top,.true.)
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top,.true.,&
+ 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_potential_acoustic,rmemory_potential_acoustic_corner,&
+ rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
+ rmemory_acoustic_dux_dx_corner,rmemory_acoustic_dux_dz_corner,deltat,&
+ PML_BOUNDARY_CONDITIONS)
endif
More information about the CIG-COMMITS
mailing list