[cig-commits] r20559 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Wed Aug 8 18:10:55 PDT 2012
Author: xie.zhinan
Date: 2012-08-08 18:10:54 -0700 (Wed, 08 Aug 2012)
New Revision: 20559
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:
fixed the some bugs for future use PML with external mesh
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-08-07 23:03:13 UTC (rev 20558)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-08-09 01:10:54 UTC (rev 20559)
@@ -55,7 +55,7 @@
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,&
- is_PML,nspec_PML,npoin_PML,ibool_PML,spec_to_PML,which_PML_elem,&
+ is_PML,nspec_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_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
@@ -127,16 +127,16 @@
integer :: ifirstelem,ilastelem
!CPML coefficients and memory variables
- integer :: nspec_PML,npoin_PML,iPML,ispec_PML
+ integer :: nspec_PML,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
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz
- 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) :: &
+ 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
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) ::PML_dux_dxl,PML_dux_dzl,&
@@ -227,16 +227,15 @@
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)) .and. &
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
!---------------------- A8 --------------------------
- A8 = - d_x_store(iPML) / (k_x_store(iPML)**2)
+ A8 = - d_x_store(i,j,ispec_PML) / (k_x_store(i,j,ispec_PML)**2)
- bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
+ bb = d_x_store(i,j,ispec_PML) / k_x_store(i,j,ispec_PML) + alpha_x_store(i,j,ispec_PML)
coef0_x = exp(-bb * deltat)
if ( abs(bb) > 1.d-3 ) then
coef1_x = (1.d0 - exp(-bb * deltat / 2.d0)) * A8 / bb
@@ -250,9 +249,9 @@
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1_x + PML_dux_dxl(i,j,ispec_PML) * coef2_x
!---------------------- A5 --------------------------
- A5 = d_x_store(iPML)
+ A5 = d_x_store(i,j,ispec_PML)
- bb = alpha_x_store(iPML)
+ bb = alpha_x_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -277,10 +276,11 @@
!------------------------------------------------------------------------------
!---------------------------- A8 ----------------------------
- A8 = (k_x_store(iPML) * d_z_store(iPML) - k_z_store(iPML) * d_x_store(iPML)) &
- / (k_x_store(iPML)**2)
+ A8 = (k_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML) &
+ - k_z_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML)) &
+ / (k_x_store(i,j,ispec_PML)**2)
- bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
+ bb = d_x_store(i,j,ispec_PML) / k_x_store(i,j,ispec_PML) + alpha_x_store(i,j,ispec_PML)
coef0_z = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -296,10 +296,11 @@
!---------------------------- A6 ----------------------------
- A6 =(k_z_store(iPML) * d_x_store(iPML) - k_x_store(iPML) * d_z_store(iPML)) &
- / (k_z_store(iPML)**2)
+ A6 =(k_z_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML) &
+ - k_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML)) &
+ / (k_z_store(i,j,ispec_PML)**2)
- bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
+ bb = d_z_store(i,j,ispec_PML) / k_z_store(i,j,ispec_PML) + alpha_z_store(i,j,ispec_PML)
coef0_z = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -322,8 +323,8 @@
!---------------------------- TOP & BOTTOM ------------------------------------
!------------------------------------------------------------------------------
!---------------------- A7 --------------------------
- A7 = d_z_store(iPML)
- bb = alpha_z_store(iPML)
+ A7 = d_z_store(i,j,ispec_PML)
+ bb = alpha_z_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -340,8 +341,8 @@
!---------------------- A6 --------------------------
- A6 = - d_z_store(iPML) / ( k_z_store(iPML) ** 2 )
- bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
+ A6 = - d_z_store(i,j,ispec_PML) / ( k_z_store(i,j,ispec_PML) ** 2 )
+ bb = d_z_store(i,j,ispec_PML) / k_z_store(i,j,ispec_PML) + alpha_z_store(i,j,ispec_PML)
coef0_x = exp(-bb * deltat)
if ( abs(bb) > 1.d-3 ) then
coef1_x = (1.d0 - exp(-bb * deltat / 2.d0)) * A6 / bb
@@ -388,7 +389,6 @@
if(is_PML(ispec) .and. PML_BOUNDARY_CONDITIONS)then
ispec_PML=spec_to_PML(ispec)
- iPML=ibool_PML(i,j,ispec)
iglob=ibool(i,j,ispec)
if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
@@ -398,9 +398,9 @@
!------------------------------------------------------------------------------
!---------------------- A3 and A4 --------------------------
- A3 = d_x_store(iPML ) * alpha_x_store(iPML) ** 2
+ A3 = d_x_store(i,j,ispec_PML) * alpha_x_store(i,j,ispec_PML) ** 2
A4 = 0.d0
- bb = alpha_x_store(iPML)
+ bb = alpha_x_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -428,7 +428,7 @@
!---------------------------- A3 ----------------------------
A3 = 1.d0
- bb = alpha_x_store(iPML)
+ bb = alpha_x_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -447,7 +447,7 @@
!---------------------------- A4 ----------------------------
A4 =1.0d0
- bb = alpha_z_store(iPML)
+ bb = alpha_z_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -471,9 +471,9 @@
!------------------------------------------------------------------------------
!---------------------- A3 and A4 ----------------------------
- A4 = d_z_store(iPML ) * alpha_z_store(iPML) ** 2
+ A4 = d_z_store(i,j,ispec_PML) * alpha_z_store(i,j,ispec_PML) ** 2
A3 = 0.d0
- bb = alpha_z_store(iPML)
+ bb = alpha_z_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -499,9 +499,9 @@
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
- A0 = - alpha_x_store( iPML ) * d_x_store(iPML)
- A1 = d_x_store(iPML)
- A2 = k_x_store(iPML)
+ A0 = - alpha_x_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML)
+ A1 = d_x_store(i,j,ispec_PML)
+ A2 = k_x_store(i,j,ispec_PML)
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
(A1*potential_dot_acoustic(iglob)+ &
@@ -514,20 +514,20 @@
!------------------------------------------------------------------------------
!-------------------------------- CORNER --------------------------------------
!------------------------------------------------------------------------------
- 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)
+ A0 = d_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML) &
+ - alpha_x_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML) &
+ - alpha_z_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)
- A1 = d_x_store(iPML) * k_z_store(iPML) + d_z_store(iPML) * k_x_store(iPML)
+ A1 = d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML) + d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)
- A2 = k_x_store(iPML) * k_z_store(iPML)
+ A2 = k_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)
- A3 = alpha_x_store(iPML) ** 2*(d_x_store(iPML) * k_z_store(iPML)+ &
- d_z_store(iPML) * k_x_store(iPML)) &
- -2.d0 * alpha_x_store(iPML)*d_x_store(iPML)*d_z_store(iPML)+ &
- (it+0.5)*deltat*alpha_x_store(iPML)**2*d_x_store(iPML)*d_z_store(iPML)
+ A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
+ d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
+ -2.d0 * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)+ &
+ (it+0.5)*deltat*alpha_x_store(i,j,ispec_PML)**2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
- A4 = -alpha_x_store(iPML) ** 2*d_x_store(iPML)*d_z_store(iPML)
+ A4 = -alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
(A1*potential_dot_acoustic(iglob)+ &
@@ -539,9 +539,9 @@
!------------------------------------------------------------------------------
!-------------------------------- TOP & BOTTOM --------------------------------
!------------------------------------------------------------------------------
- A0 = - alpha_z_store( iPML ) * d_z_store(iPML)
- A1 = d_z_store(iPML)
- A2 = k_z_store(iPML)
+ A0 = - alpha_z_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML)
+ A1 = d_z_store(i,j,ispec_PML)
+ A2 = k_z_store(i,j,ispec_PML)
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
(A1*potential_dot_acoustic(iglob)+ &
@@ -573,7 +573,6 @@
if(is_PML(ispec) .and. PML_BOUNDARY_CONDITIONS)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)) .and. &
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-08-07 23:03:13 UTC (rev 20558)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-08-09 01:10:54 UTC (rev 20559)
@@ -63,7 +63,7 @@
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,nadj_rec_local, &
- is_PML,nspec_PML,npoin_PML,ibool_PML,spec_to_PML,which_PML_elem, &
+ is_PML,nspec_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_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
PML_BOUNDARY_CONDITIONS)
@@ -214,17 +214,17 @@
integer :: ifirstelem,ilastelem
!CPML coefficients and memory variables
- integer :: nspec_PML,npoin_PML,iPML,ispec_PML
+ integer :: nspec_PML,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
logical :: PML_BOUNDARY_CONDITIONS
real(kind=CUSTOM_REAL), dimension(2,3,NGLLX,NGLLZ,nspec_PML) :: rmemory_displ_elastic
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz
- 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) :: &
+ 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
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
@@ -408,7 +408,6 @@
if(is_PML(ispec) .and. PML_BOUNDARY_CONDITIONS) then
ispec_PML=spec_to_PML(ispec)
- iPML=ibool_PML(i,j,ispec)
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
@@ -416,8 +415,8 @@
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
!---------------------- A8--------------------------
- A8 = - d_x_store(iPML) / (k_x_store(iPML) ** 2)
- bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
+ A8 = - d_x_store(i,j,ispec_PML) / (k_x_store(i,j,ispec_PML) ** 2)
+ bb = d_x_store(i,j,ispec_PML) / k_x_store(i,j,ispec_PML) + alpha_x_store(i,j,ispec_PML)
coef0_x = exp(-bb * deltat)
if ( abs(bb) > 1.d-3 ) then
coef1_x = (1.d0 - exp(-bb * deltat / 2.d0)) * A8 / bb
@@ -435,8 +434,8 @@
+ PML_duz_dxl_new(i,j,ispec_PML) * coef1_x + PML_duz_dxl(i,j,ispec_PML) * coef2_x
!---------------------- A5--------------------------
- A5 = d_x_store(iPML)
- bb = alpha_x_store(iPML)
+ A5 = d_x_store(i,j,ispec_PML)
+ bb = alpha_x_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -467,10 +466,11 @@
elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
(which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
- A8 = (k_x_store(iPML) * d_z_store(iPML) - k_z_store(iPML) * d_x_store(iPML)) &
- / (k_x_store(iPML)**2)
+ A8 = (k_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML) &
+ - k_z_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML)) &
+ / (k_x_store(i,j,ispec_PML)**2)
- bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
+ bb = d_x_store(i,j,ispec_PML) / k_x_store(i,j,ispec_PML) + alpha_x_store(i,j,ispec_PML)
coef0_z = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -490,10 +490,11 @@
+ PML_duz_dxl_new(i,j,ispec_PML) * coef1_z + PML_duz_dxl(i,j,ispec_PML) * coef2_z
!---------------------------- A6 ----------------------------
- A6 =(k_z_store(iPML) * d_x_store(iPML) - k_x_store(iPML) * d_z_store(iPML)) &
- / (k_z_store(iPML)**2)
+ A6 =(k_z_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML) &
+ - k_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML)) &
+ / (k_z_store(i,j,ispec_PML)**2)
- bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
+ bb = d_z_store(i,j,ispec_PML) / k_z_store(i,j,ispec_PML) + alpha_z_store(i,j,ispec_PML)
coef0_z = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -523,8 +524,8 @@
!------------------------------------------------------------------------------
!---------------------- A7 --------------------------
- A7 = d_z_store(iPML)
- bb = alpha_z_store(iPML)
+ A7 = d_z_store(i,j,ispec_PML)
+ bb = alpha_z_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -545,8 +546,8 @@
+ PML_duz_dxl_new(i,j,ispec_PML) * coef1_x + PML_duz_dxl(i,j,ispec_PML) * coef2_x
!---------------------- A6 --------------------------
- A6 = - d_z_store(iPML) / ( k_z_store(iPML) ** 2 )
- bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
+ A6 = - d_z_store(i,j,ispec_PML) / ( k_z_store(i,j,ispec_PML) ** 2 )
+ bb = d_z_store(i,j,ispec_PML) / k_z_store(i,j,ispec_PML) + alpha_z_store(i,j,ispec_PML)
coef0_x = exp(-bb * deltat)
if ( abs(bb) > 1.d-3 ) then
coef1_x = (1.d0 - exp(-bb * deltat / 2.d0)) * A6 / bb
@@ -649,7 +650,6 @@
if(PML_BOUNDARY_CONDITIONS .and. 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)
@@ -755,7 +755,6 @@
else
rhol = density(1,kmato(ispec))
endif
- iPML=ibool_PML(i,j,ispec)
iglob=ibool(i,j,ispec)
!------------------------------------------------------------------------------
@@ -765,9 +764,9 @@
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
!---------------------- A3 and A4 --------------------------
- A3 = d_x_store(iPML ) * alpha_x_store(iPML) ** 2
+ A3 = d_x_store(i,j,ispec_PML) * alpha_x_store(i,j,ispec_PML) ** 2
A4 = 0.d0
- bb = alpha_x_store(iPML)
+ bb = alpha_x_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -797,7 +796,7 @@
A3 = 1.d0
- bb = alpha_x_store(iPML)
+ bb = alpha_x_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -818,7 +817,7 @@
!---------------------------- A4 ----------------------------
A4 =1.0d0
- bb = alpha_z_store(iPML)
+ bb = alpha_z_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -842,9 +841,9 @@
!------------------------------------------------------------------------------
!---------------------- A3 and A4 ----------------------------
- A4 = d_z_store(iPML ) * alpha_z_store(iPML) ** 2
+ A4 = d_z_store(i,j,ispec_PML) * alpha_z_store(i,j,ispec_PML) ** 2
A3 = 0.d0
- bb = alpha_z_store(iPML)
+ bb = alpha_z_store(i,j,ispec_PML)
coef0_x = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -870,9 +869,9 @@
if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
- A0 = - alpha_x_store( iPML ) * d_x_store(iPML)
- A1 = d_x_store(iPML)
- A2 = k_x_store(iPML)
+ A0 = - alpha_x_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML)
+ A1 = d_x_store(i,j,ispec_PML)
+ A2 = k_x_store(i,j,ispec_PML)
accel_elastic_PML(1,i,j,ispec_PML) = wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( A1 * veloc_elastic(1,iglob)+ &
@@ -887,19 +886,19 @@
elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
(which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
- 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)
+ A0 = d_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML) &
+ - alpha_x_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML) &
+ - alpha_z_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)
- A1 = d_x_store(iPML) * k_z_store(iPML) + d_z_store(iPML) * k_x_store(iPML)
+ A1 = d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML) + d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)
- A2 = k_x_store(iPML) * k_z_store(iPML)
+ A2 = k_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)
- A3 = alpha_x_store(iPML) ** 2*(d_x_store(iPML) * k_z_store(iPML)+ &
- d_z_store(iPML) * k_x_store(iPML)) &
- -2.d0 * alpha_x_store(iPML)*d_x_store(iPML)*d_z_store(iPML)+ &
- (it+0.5)*deltat*alpha_x_store(iPML)**2*d_x_store(iPML)*d_z_store(iPML)
- A4 = -alpha_x_store(iPML) ** 2*d_x_store(iPML)*d_z_store(iPML)
+ A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
+ d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
+ -2.d0 * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)+ &
+ (it+0.5)*deltat*alpha_x_store(i,j,ispec_PML)**2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
+ A4 = -alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
accel_elastic_PML(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( A1 *veloc_elastic(1,iglob)+ &
@@ -913,9 +912,9 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
else
- A0 = - alpha_z_store( iPML ) * d_z_store(iPML)
- A1 = d_z_store(iPML)
- A2 = k_z_store(iPML)
+ A0 = - alpha_z_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML)
+ A1 = d_z_store(i,j,ispec_PML)
+ A2 = k_z_store(i,j,ispec_PML)
accel_elastic_PML(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( A1 * veloc_elastic(1,iglob)+ &
@@ -964,7 +963,6 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(is_PML(ispec) .and. PML_BOUNDARY_CONDITIONS)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)) .and. &
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
accel_elastic(1,iglob) = accel_elastic(1,iglob) - accel_elastic_PML(1,i,j,ispec_PML)
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-08-07 23:03:13 UTC (rev 20558)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-08-09 01:10:54 UTC (rev 20559)
@@ -59,8 +59,9 @@
#ifdef USE_GUENNEAU
,coord &
#endif
- ,K_x_store,K_z_store,npoin_PML,ibool_PML,is_PML,&
- d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,which_PML_elem)
+ ,K_x_store,K_z_store,is_PML,&
+ d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,which_PML_elem,&
+ nspec_PML,spec_to_PML)
! builds the global mass matrix
@@ -123,9 +124,10 @@
double precision :: rhol_s,rhol_f,rhol_bar,phil,tortl
!!!!!!!!!!!!! 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
- integer, dimension(NGLLX,NGLLZ,nspec) :: ibool_PML
+ integer :: nspec_PML,ispec_PML
+ integer, dimension(nspec) :: spec_to_PML
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
+ K_x_store,K_z_store,d_x_store,d_z_store
logical, dimension(nspec) :: is_PML
logical :: PML_BOUNDARY_CONDITIONS
logical, dimension(4,nspec) :: which_PML_elem
@@ -182,23 +184,24 @@
! for elastic medium
if (is_PML(ispec) .and. PML_BOUNDARY_CONDITIONS) then
- iPML=ibool_PML(i,j,ispec)
+ ispec_PML=spec_to_PML(ispec)
if ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
- + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(iPML)&
- + d_x_store(iPML) * deltat / 2.d0)
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML)&
+ + d_x_store(i,j,ispec_PML) * deltat / 2.d0)
rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
(which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
- + wxgll(i)*wzgll(j)*rhol*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)
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML) * K_z_store(i,j,ispec_PML)&
+ + (d_x_store(i,j,ispec_PML)*k_z_store(i,j,ispec_PML)+&
+ d_z_store(i,j,ispec_PML)*k_x_store(i,j,ispec_PML)) * deltat / 2.d0)
rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
else
rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
- + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_z_store(iPML)&
- + d_z_store(iPML)* deltat / 2.d0)
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_z_store(i,j,ispec_PML)&
+ + d_z_store(i,j,ispec_PML)* deltat / 2.d0)
rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
endif
@@ -230,22 +233,22 @@
! + wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / kappal
if (PML_BOUNDARY_CONDITIONS .and. is_PML(ispec)) then
- iPML=ibool_PML(i,j,ispec)
-
+ ispec_PML=spec_to_PML(ispec)
if ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
- + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(iPML)&
- + d_x_store(iPML) * deltat / 2.d0)
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML)&
+ + d_x_store(i,j,ispec_PML) * deltat / 2.d0)
elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
(which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
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)
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML) * K_z_store(i,j,ispec_PML)&
+ + (d_x_store(i,j,ispec_PML)*k_z_store(i,j,ispec_PML)&
+ +d_z_store(i,j,ispec_PML)*k_x_store(i,j,ispec_PML)) * deltat / 2.d0)
else
rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
- + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_z_store(iPML)&
- + d_z_store(iPML)* deltat / 2.d0)
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_z_store(i,j,ispec_PML)&
+ + d_z_store(i,j,ispec_PML)* deltat / 2.d0)
endif
else
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-08-07 23:03:13 UTC (rev 20558)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-08-09 01:10:54 UTC (rev 20559)
@@ -42,8 +42,8 @@
!========================================================================
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,NELEM_PML_THICKNESS)
+ nspec_PML,is_PML,which_PML_elem,which_PML_poin,spec_to_PML, &
+ icorner_iglob,NELEM_PML_THICKNESS)
implicit none
@@ -53,7 +53,7 @@
include 'mpif.h'
#endif
- integer :: nspec,nglob,nelemabs,nspec_PML,npoin_PML,NELEM_PML_THICKNESS
+ integer :: nspec,nglob,nelemabs,nspec_PML,NELEM_PML_THICKNESS
logical :: anyabs
integer :: ibound,ispecabs,ncorner,ispec,iglob
@@ -66,8 +66,6 @@
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
!!!detection of PML elements
@@ -158,7 +156,6 @@
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
@@ -168,40 +165,7 @@
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)
-
write(IOUT,*) "number of PML spectral elements :", nspec_PML
- write(IOUT,*) "number of PML spectral points :", npoin_PML
end subroutine pml_init
@@ -210,9 +174,10 @@
!
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)
+ which_PML_elem,kmato,density,poroelastcoef,numat,f0_temp,&
+ myrank,&
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store,&
+ nspec_PML,spec_to_PML)
implicit none
@@ -222,16 +187,16 @@
include "mpif.h"
#endif
- integer nspec, npoin, i, j,numat, ispec,iglob,npoin_PML,iPML
+ integer nspec, npoin, i, j,numat, ispec,iglob,nspec_PML
double precision :: f0_temp
logical, dimension(nspec) :: is_PML
logical, dimension(4,nspec) :: which_PML_elem
- real(kind=CUSTOM_REAL), dimension(npoin_PML) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
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
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
double precision, dimension(2,numat) :: density
double precision, dimension(4,3,numat) :: poroelastcoef
integer, dimension(nspec) :: kmato
@@ -275,11 +240,13 @@
double precision :: Rcoef
- double precision :: factorx, factorz
-
! material properties of the elastic medium
double precision :: lambdalplus2mul_relaxed
+!DK,DK
+ integer, dimension(nspec) :: spec_to_PML
+ integer :: ispec_PML
+
! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
@@ -441,14 +408,18 @@
do ispec = 1,nspec
-
+ ispec_PML=spec_to_PML(ispec)
if (is_PML(ispec)) then
-
do j=1,NGLLZ
do i=1,NGLLX
+ K_x_store(i,j,ispec_PML) = 0.0
+ K_z_store(i,j,ispec_PML) = 0.0
+ d_x_store(i,j,ispec_PML) = 0.0
+ d_z_store(i,j,ispec_PML) = 0.0
+ alpha_x_store(i,j,ispec_PML) = 0.0
+ alpha_z_store(i,j,ispec_PML) = 0.0
+
iglob=ibool(i,j,ispec)
- iPML=ibool_PML(i,j,ispec)
- if(iPML < 1) stop 'error: iPML < 1 in a PML element'
! abscissa of current grid point along the damping profile
xval = coord(1,iglob)
zval = coord(2,iglob)
@@ -479,12 +450,6 @@
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
@@ -507,12 +472,6 @@
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
@@ -535,12 +494,6 @@
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
@@ -563,20 +516,14 @@
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
- K_x_store(iPML) = K_x
- K_z_store(iPML) = K_z
- d_x_store(iPML) = d_x
- d_z_store(iPML) = d_z
- alpha_x_store(iPML) = alpha_x
- alpha_z_store(iPML) = alpha_z
+ K_x_store(i,j,ispec_PML) = K_x
+ K_z_store(i,j,ispec_PML) = K_z
+ d_x_store(i,j,ispec_PML) = d_x
+ d_z_store(i,j,ispec_PML) = d_z
+ alpha_x_store(i,j,ispec_PML) = alpha_x
+ alpha_z_store(i,j,ispec_PML) = alpha_z
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-08-07 23:03:13 UTC (rev 20558)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-08-09 01:10:54 UTC (rev 20559)
@@ -989,14 +989,14 @@
!PML parameters
logical, dimension(:), allocatable :: is_PML
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
+
+ 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
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
@@ -1006,7 +1006,7 @@
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz
logical :: anyabs_glob
- integer :: nspec_PML, npoin_PML
+ integer :: nspec_PML
!***********************************************************************
!
@@ -2817,43 +2817,43 @@
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,NELEM_PML_THICKNESS)
+ nspec_PML,is_PML,which_PML_elem,which_PML_poin,spec_to_PML, &
+ icorner_iglob,NELEM_PML_THICKNESS)
deallocate(icorner_iglob)
deallocate(which_PML_poin)
- if (npoin_PML==0) npoin_PML=1
+ if (nspec_PML==0) nspec_PML=1
!!!!!!!!!!!!! DK DK added this
- if (npoin_PML > 0) then
+ if (nspec_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))
+ allocate(K_x_store(NGLLX,NGLLZ,nspec_PML))
+ allocate(K_z_store(NGLLX,NGLLZ,nspec_PML))
+ allocate(d_x_store(NGLLX,NGLLZ,nspec_PML))
+ allocate(d_z_store(NGLLX,NGLLZ,nspec_PML))
+ allocate(alpha_x_store(NGLLX,NGLLZ,nspec_PML))
+ allocate(alpha_z_store(NGLLX,NGLLZ,nspec_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
+ 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)
+ which_PML_elem,kmato,density,poroelastcoef,numat,f0(1),&
+ myrank,&
+ K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store,&
+ nspec_PML,spec_to_PML)
endif
@@ -2923,16 +2923,15 @@
allocate(rmemory_acoustic_dux_dz(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))
+ allocate(K_x_store(1,1,1))
+ allocate(K_z_store(1,1,1))
+ allocate(d_x_store(1,1,1))
+ allocate(d_z_store(1,1,1))
+ allocate(alpha_x_store(1,1,1))
+ allocate(alpha_z_store(1,1,1))
end if ! PML_BOUNDARY_CONDITIONS
@@ -2959,8 +2958,9 @@
#ifdef USE_GUENNEAU
,coord &
#endif
- ,K_x_store,K_z_store,npoin_PML,ibool_PML,is_PML,&
- d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,which_PML_elem)
+ ,K_x_store,K_z_store,is_PML,&
+ d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,which_PML_elem,&
+ nspec_PML,spec_to_PML)
#ifdef USE_MPI
if ( nproc > 1 ) then
@@ -4879,7 +4879,7 @@
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.,&
- is_PML,nspec_PML,npoin_PML,ibool_PML,spec_to_PML,which_PML_elem,&
+ is_PML,nspec_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_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
@@ -4898,7 +4898,7 @@
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.,&
- is_PML,nspec_PML,npoin_PML,ibool_PML,spec_to_PML,which_PML_elem,&
+ is_PML,nspec_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_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
@@ -5373,7 +5373,7 @@
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,max(1,nadj_rec_local), &
- is_PML,nspec_PML,npoin_PML,ibool_PML,spec_to_PML,which_PML_elem, &
+ is_PML,nspec_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_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz,&
PML_BOUNDARY_CONDITIONS)
More information about the CIG-COMMITS
mailing list