[cig-commits] r20535 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon Jul 23 05:25:31 PDT 2012
Author: xie.zhinan
Date: 2012-07-23 05:25:30 -0700 (Mon, 23 Jul 2012)
New Revision: 20535
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/pml_init.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
Simplfied the code by setting alpha_x=alpha_z=const where alpha_x or alpha_z is nonzero
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-07-22 15:46:42 UTC (rev 20534)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-07-23 12:25:30 UTC (rev 20535)
@@ -134,7 +134,7 @@
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(2,NGLLX,NGLLZ,nspec_PML) :: &
+ 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
@@ -233,9 +233,8 @@
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
- !---------------------- A8 and A6 --------------------------
- A8 = - d_x_store(iPML) / (k_x_store(iPML) ** 2)
- A6 = 0.d0
+ !---------------------- A8 --------------------------
+ A8 = - d_x_store(iPML) / (k_x_store(iPML)**2)
bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
coef0_x = exp(-bb * deltat)
@@ -247,13 +246,11 @@
coef2_x = deltat * A8
end if
- 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) &
+ rmemory_acoustic_dux_dx(i,j,ispec_PML) = coef0_x*rmemory_acoustic_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1_x + PML_dux_dxl(i,j,ispec_PML) * coef2_x
- !---------------------- A5 and A7 --------------------------
+ !---------------------- A5 --------------------------
A5 = d_x_store(iPML)
- A7 = 0.d0
bb = alpha_x_store(iPML)
coef0_x = exp(- bb * deltat)
@@ -267,76 +264,22 @@
coef2_x = deltat * A5 ! Rene Matzen
end if
- rmemory_acoustic_dux_dz(1,i,j,ispec_PML) = coef0_x * rmemory_acoustic_dux_dz(1,i,j,ispec_PML) &
+ rmemory_acoustic_dux_dz(i,j,ispec_PML) = coef0_x * rmemory_acoustic_dux_dz(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)
+ dux_dxl = PML_dux_dxl(i,j,ispec_PML) + rmemory_acoustic_dux_dx(i,j,ispec_PML)
+ dux_dzl = PML_dux_dzl(i,j,ispec_PML) + rmemory_acoustic_dux_dz(i,j,ispec_PML)
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
!------------------------------------------------------------------------------
!---------------------------- CORNER ------------------------------------------
!------------------------------------------------------------------------------
- !----------------------- 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
-
- 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
-
-
!---------------------------- 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
+ A8 = (k_x_store(iPML) * d_z_store(iPML) - k_z_store(iPML) * d_x_store(iPML)) &
+ / (k_x_store(iPML)**2)
+
bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
coef0_z = exp(- bb * deltat)
@@ -348,66 +291,14 @@
coef2_z = deltat * A8 ! Rene Matzen
end if
- rmemory_acoustic_dux_dx(2,i,j,ispec_PML) = coef0_z*rmemory_acoustic_dux_dx(2,i,j,ispec_PML) &
+ rmemory_acoustic_dux_dx(i,j,ispec_PML) = coef0_z*rmemory_acoustic_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1_z + PML_dux_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
-
- 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
-
-
!---------------------------- 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
+ A6 =(k_z_store(iPML) * d_x_store(iPML) - k_x_store(iPML) * d_z_store(iPML)) &
+ / (k_z_store(iPML)**2)
+
bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
coef0_z = exp(- bb * deltat)
@@ -419,22 +310,19 @@
coef2_z = deltat * A6 ! Rene Matzen
end if
- rmemory_acoustic_dux_dz(2,i,j,ispec_PML) = coef0_z * rmemory_acoustic_dux_dz(2,i,j,ispec_PML) &
+ rmemory_acoustic_dux_dz(i,j,ispec_PML) = coef0_z * rmemory_acoustic_dux_dz(i,j,ispec_PML) &
+ PML_dux_dzl_new(i,j,ispec_PML) *coef1_z + PML_dux_dzl(i,j,ispec_PML) * coef2_z
- 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)
+ dux_dxl = PML_dux_dxl(i,j,ispec_PML) + rmemory_acoustic_dux_dx(i,j,ispec_PML)
+ dux_dzl = PML_dux_dzl(i,j,ispec_PML) + rmemory_acoustic_dux_dz(i,j,ispec_PML)
else
!------------------------------------------------------------------------------
!---------------------------- TOP & BOTTOM ------------------------------------
!------------------------------------------------------------------------------
- !---------------------- A5 and A7 --------------------------
+ !---------------------- A7 --------------------------
A7 = d_z_store(iPML)
- A5 = 0.d0
bb = alpha_z_store(iPML)
coef0_x = exp(- bb * deltat)
@@ -447,14 +335,12 @@
coef2_x = deltat * A7 ! Rene Matzen
end if
- rmemory_acoustic_dux_dx(1,i,j,ispec_PML) = coef0_x*rmemory_acoustic_dux_dx(1,i,j,ispec_PML) &
+ rmemory_acoustic_dux_dx(i,j,ispec_PML) = coef0_x*rmemory_acoustic_dux_dx(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
- !---------------------- A8 and A6 --------------------------
+ !---------------------- 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
@@ -465,12 +351,11 @@
coef2_x = deltat * A6
end if
- 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) &
+ rmemory_acoustic_dux_dz(i,j,ispec_PML) = coef0_x * rmemory_acoustic_dux_dz(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)
+ dux_dxl = dux_dxl + rmemory_acoustic_dux_dx(i,j,ispec_PML)
+ dux_dzl = dux_dzl + rmemory_acoustic_dux_dz(i,j,ispec_PML)
endif
endif
@@ -541,22 +426,8 @@
!------------------------------------------------------------------------------
!---------------------------- 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 &
- ) then
- 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) ) &
- / ( alpha_x_store(iPML) - alpha_z_store(iPML) )
- end if
+ A3 = 1.d0
+
bb = alpha_x_store(iPML)
coef0_x = exp(- bb * deltat)
@@ -574,22 +445,8 @@
+ potential_acoustic(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 &
- ) then
- 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) ) &
- / ( alpha_z_store(iPML) - alpha_x_store(iPML) )
- end if
+ A4 =1.0d0
+
bb = alpha_z_store(iPML)
coef0_x = exp(- bb * deltat)
@@ -603,8 +460,8 @@
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
+ + (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob))*(it+0.5)*deltat * coef1_x &
+ + potential_acoustic(iglob) *(it-0.5)*deltat * coef2_x
else
@@ -665,9 +522,16 @@
A2 = k_x_store(iPML) * k_z_store(iPML)
+ 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)
+
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
(A1*potential_dot_acoustic(iglob)+ &
- rmemory_potential_acoustic(1,i,j,ispec_PML)+rmemory_potential_acoustic(2,i,j,ispec_PML)&
+ A3*rmemory_potential_acoustic(1,i,j,ispec_PML)+A4*rmemory_potential_acoustic(2,i,j,ispec_PML)&
+A0*potential_acoustic(iglob))
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-07-22 15:46:42 UTC (rev 20534)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-07-23 12:25:30 UTC (rev 20535)
@@ -222,7 +222,7 @@
logical :: PML_BOUNDARY_CONDITIONS
real(kind=CUSTOM_REAL), dimension(2,3,NGLLX,NGLLZ,nspec_PML) :: rmemory_displ_elastic
- real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLZ,nspec_PML) :: &
+ 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(3,NGLLX,NGLLZ,nspec_PML) :: accel_elastic_PML
@@ -415,10 +415,8 @@
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
- !---------------------- A8 and A6 --------------------------
+ !---------------------- A8--------------------------
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
@@ -428,21 +426,16 @@
coef1_x = deltat / 2.0d0 * A8
coef2_x = deltat * A8
end if
-
- 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) &
+ rmemory_dux_dx(i,j,ispec_PML) = coef0_x*rmemory_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1_x + PML_dux_dxl(i,j,ispec_PML) * coef2_x
- 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) &
+ rmemory_duz_dx(i,j,ispec_PML) = coef0_x*rmemory_duz_dx(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--------------------------
A5 = d_x_store(iPML)
- A7 = 0.d0
-
bb = alpha_x_store(iPML)
coef0_x = exp(- bb * deltat)
@@ -456,19 +449,17 @@
end if
!! 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) &
+ rmemory_dux_dz(i,j,ispec_PML) = coef0_x * rmemory_dux_dz(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
!! 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) &
+ rmemory_duz_dz(i,j,ispec_PML) = coef0_x * rmemory_duz_dz(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)
+ dux_dxl = PML_dux_dxl(i,j,ispec_PML) + rmemory_dux_dx(i,j,ispec_PML)
+ duz_dxl = PML_duz_dxl(i,j,ispec_PML) + rmemory_duz_dx(i,j,ispec_PML)
+ dux_dzl = PML_dux_dzl(i,j,ispec_PML) + rmemory_dux_dz(i,j,ispec_PML)
+ duz_dzl = PML_duz_dzl(i,j,ispec_PML) + rmemory_duz_dz(i,j,ispec_PML)
!------------------------------------------------------------------------------
!---------------------------- CORNER ------------------------------------------
@@ -476,66 +467,9 @@
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
- !----------------------- 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)
+ A8 = (k_x_store(iPML) * d_z_store(iPML) - k_z_store(iPML) * d_x_store(iPML)) &
+ / (k_x_store(iPML)**2)
- 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(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
-
- !! 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
-
- !---------------------------- 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)
@@ -548,73 +482,17 @@
end if
!! DK DK new from Wang eq (21)
- rmemory_dux_dx(2,i,j,ispec_PML) = coef0_z*rmemory_dux_dx(2,i,j,ispec_PML) &
+ rmemory_dux_dx(i,j,ispec_PML) = coef0_z*rmemory_dux_dx(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(2,i,j,ispec_PML) = coef0_z*rmemory_duz_dx(2,i,j,ispec_PML) &
+ rmemory_duz_dx(i,j,ispec_PML) = coef0_z*rmemory_duz_dx(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(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
-
- !! 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
-
!---------------------------- 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
+ A6 =(k_z_store(iPML) * d_x_store(iPML) - k_x_store(iPML) * d_z_store(iPML)) &
+ / (k_z_store(iPML)**2)
+
bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
coef0_z = exp(- bb * deltat)
@@ -627,30 +505,25 @@
end if
!! DK DK new from Wang eq (21)
- rmemory_dux_dz(2,i,j,ispec_PML) = coef0_z * rmemory_dux_dz(2,i,j,ispec_PML) &
+ rmemory_dux_dz(i,j,ispec_PML) = coef0_z * rmemory_dux_dz(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(2,i,j,ispec_PML) = coef0_z * rmemory_duz_dz(2,i,j,ispec_PML) &
+ rmemory_duz_dz(i,j,ispec_PML) = coef0_z * rmemory_duz_dz(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(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)
+ dux_dxl = PML_dux_dxl(i,j,ispec_PML) + rmemory_dux_dx(i,j,ispec_PML)
+ duz_dxl = PML_duz_dxl(i,j,ispec_PML) + rmemory_duz_dx(i,j,ispec_PML)
+ dux_dzl = PML_dux_dzl(i,j,ispec_PML) + rmemory_dux_dz(i,j,ispec_PML)
+ duz_dzl = PML_duz_dzl(i,j,ispec_PML) + rmemory_duz_dz(i,j,ispec_PML)
else
!------------------------------------------------------------------------------
!---------------------------- TOP & BOTTOM ------------------------------------
!------------------------------------------------------------------------------
- !---------------------- A5 and A7 --------------------------
+ !---------------------- A7 --------------------------
A7 = d_z_store(iPML)
- A5 = 0.d0
bb = alpha_z_store(iPML)
coef0_x = exp(- bb * deltat)
@@ -664,18 +537,15 @@
end if
!! 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) &
+ rmemory_dux_dx(i,j,ispec_PML) = coef0_x*rmemory_dux_dx(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
!! 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) &
+ rmemory_duz_dx(i,j,ispec_PML) = coef0_x*rmemory_duz_dx(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 --------------------------
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
@@ -686,20 +556,18 @@
coef2_x = deltat * A6
end if
- 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) &
+ rmemory_dux_dz(i,j,ispec_PML) = coef0_x * rmemory_dux_dz(i,j,ispec_PML) &
+ PML_dux_dzl_new(i,j,ispec_PML) *coef1_x + PML_dux_dzl(i,j,ispec_PML) * coef2_x
- 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) &
+ rmemory_duz_dz(i,j,ispec_PML) = coef0_x * rmemory_duz_dz(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)
+ dux_dxl = dux_dxl + rmemory_dux_dx(i,j,ispec_PML)
+ duz_dxl = duz_dxl + rmemory_duz_dx(i,j,ispec_PML)
+ dux_dzl = dux_dzl + rmemory_dux_dz(i,j,ispec_PML)
+ duz_dzl = duz_dzl + rmemory_duz_dz(i,j,ispec_PML)
endif
endif ! PML_BOUNDARY_CONDITIONS
@@ -926,22 +794,9 @@
(which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
!---------------------------- 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 &
- ) then
- 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) ) &
- / ( alpha_x_store(iPML) - alpha_z_store(iPML) )
- end if
+
+ A3 = 1.d0
+
bb = alpha_x_store(iPML)
coef0_x = exp(- bb * deltat)
@@ -952,6 +807,7 @@
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
@@ -959,24 +815,9 @@
coef0_x * rmemory_displ_elastic(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 &
- ) then
- 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) ) &
- / ( alpha_z_store(iPML) - alpha_x_store(iPML) )
- end if
+ A4 =1.0d0
+
bb = alpha_z_store(iPML)
coef0_x = exp(- bb * deltat)
@@ -988,12 +829,12 @@
coef2_x = deltat * A4 ! Rene Matzen
end if
- 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(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
+ 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)*(it+0.5)*deltat * coef1_x &
+ + displ_elastic(1,iglob)*(it-0.5)*deltat * coef2_x
+ 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)*(it+0.5)*deltat * coef1_x &
+ + displ_elastic(3,iglob)*(it-0.5)*deltat * coef2_x
else
!------------------------------------------------------------------------------
@@ -1054,13 +895,19 @@
A2 = k_x_store(iPML) * k_z_store(iPML)
+ 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)
+
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)&
+ A3*rmemory_displ_elastic(1,1,i,j,ispec_PML)+A4*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)&
+ A3*rmemory_displ_elastic(1,3,i,j,ispec_PML)+A4*rmemory_displ_elastic(2,3,i,j,ispec_PML)&
+ A0 * displ_elastic(3,iglob) )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-07-22 15:46:42 UTC (rev 20534)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-07-23 12:25:30 UTC (rev 20535)
@@ -459,11 +459,20 @@
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) &
- + ALPHA_MAX_PML / 2.d0
+!DK DK we keep the equation to define an nonconstant alpha_z in case user needed,
+!However for users who want to use nonconstant alpha_z, you also need to change
+!routines for CMPL computation. For example in compute_forces_viscoelastic.f90
+! alpha_z = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
+! + ALPHA_MAX_PML / 2.d0
+
+!DK DK Here we set alpha_z=alpha_x=const where alpha_z or alpha_x is nonzero
+
+ alpha_z = ALPHA_MAX_PML / 2.d0
+
K_z = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
else
d_z = 0.d0
@@ -484,10 +493,14 @@
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) &
- + ALPHA_MAX_PML / 2.d0
+! alpha_z = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
+! + ALPHA_MAX_PML / 2.d0
+
+ alpha_z = ALPHA_MAX_PML / 2.d0
+
K_z = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
else
d_z = 0.d0
@@ -508,10 +521,14 @@
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
- alpha_x = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
- + ALPHA_MAX_PML / 2.d0
+! alpha_x = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
+! + ALPHA_MAX_PML / 2.d0
+
+ alpha_x = ALPHA_MAX_PML / 2.d0
+
K_x = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
else
d_x = 0.d0
@@ -532,10 +549,14 @@
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
- alpha_x = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
- + ALPHA_MAX_PML / 2.d0
+! alpha_x = ALPHA_MAX_PML * (1.d0 - abscissa_normalized) &
+! + ALPHA_MAX_PML / 2.d0
+
+ alpha_x = ALPHA_MAX_PML / 2.d0
+
K_x = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
else
d_x = 0.d0
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-07-22 15:46:42 UTC (rev 20534)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-07-23 12:25:30 UTC (rev 20535)
@@ -992,7 +992,7 @@
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 :: &
+ 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
@@ -1001,7 +1001,8 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rmemory_potential_acoustic,&
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rmemory_potential_acoustic
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz
logical :: anyabs_glob
@@ -2861,28 +2862,28 @@
allocate(rmemory_displ_elastic(2,3,NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_displ_elastic'
- allocate(rmemory_dux_dx(2,NGLLX,NGLLZ,nspec_PML),stat=ier)
+ allocate(rmemory_dux_dx(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_dux_dx'
- allocate(rmemory_dux_dz(2,NGLLX,NGLLZ,nspec_PML),stat=ier)
+ allocate(rmemory_dux_dz(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_dux_dz'
- allocate(rmemory_duz_dx(2,NGLLX,NGLLZ,nspec_PML),stat=ier)
+ allocate(rmemory_duz_dx(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_duz_dx'
- allocate(rmemory_duz_dz(2,NGLLX,NGLLZ,nspec_PML),stat=ier)
+ allocate(rmemory_duz_dz(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_duz_dz'
rmemory_displ_elastic(:,:,:,:,:) = ZERO
- rmemory_dux_dx(:,:,:,:) = ZERO
- rmemory_dux_dz(:,:,:,:) = ZERO
- rmemory_duz_dx(:,:,:,:) = ZERO
- rmemory_duz_dz(:,:,:,:) = ZERO
+ rmemory_dux_dx(:,:,:) = ZERO
+ rmemory_dux_dz(:,:,:) = ZERO
+ rmemory_duz_dx(:,:,:) = ZERO
+ rmemory_duz_dz(:,:,:) = ZERO
else
allocate(rmemory_displ_elastic(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(1,1,1))
+ allocate(rmemory_dux_dz(1,1,1))
+ allocate(rmemory_duz_dx(1,1,1))
+ allocate(rmemory_duz_dz(1,1,1))
end if
@@ -2890,9 +2891,9 @@
allocate(rmemory_potential_acoustic(2,NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_potential_acoustic'
- allocate(rmemory_acoustic_dux_dx(2,NGLLX,NGLLZ,nspec_PML),stat=ier)
+ allocate(rmemory_acoustic_dux_dx(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_acoustic_dux_dx'
- allocate(rmemory_acoustic_dux_dz(2,NGLLX,NGLLZ,nspec_PML),stat=ier)
+ allocate(rmemory_acoustic_dux_dz(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_acoustic_dux_dz'
rmemory_potential_acoustic = ZERO
@@ -2902,8 +2903,8 @@
else
allocate(rmemory_potential_acoustic(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(1,1,1))
+ allocate(rmemory_acoustic_dux_dz(1,1,1))
end if
@@ -2912,14 +2913,14 @@
!!!!!!!!!!!!! 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(1,1,1))
+ allocate(rmemory_dux_dz(1,1,1))
+ allocate(rmemory_duz_dx(1,1,1))
+ allocate(rmemory_duz_dz(1,1,1))
allocate(rmemory_displ_elastic(1,1,1,1,1))
allocate(rmemory_potential_acoustic(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(1,1,1))
+ allocate(rmemory_acoustic_dux_dz(1,1,1))
allocate(is_PML(1))
allocate(ibool_PML(1,1,1))
More information about the CIG-COMMITS
mailing list