[cig-commits] r20522 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Fri Jul 13 01:58:42 PDT 2012
Author: xie.zhinan
Date: 2012-07-13 01:58:41 -0700 (Fri, 13 Jul 2012)
New Revision: 20522
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
clean the code compute_forces_acoustic.f90, add some stop flag to prevent some errors come from allocate dimension
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-07-11 21:26:25 UTC (rev 20521)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-07-13 08:58:41 UTC (rev 20522)
@@ -141,14 +141,22 @@
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,&
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
+ real(kind=CUSTOM_REAL) :: coef0_x, coef1_x, coef2_x, coef0_z, coef1_z, coef2_z,bb,deltat
+ real(kind=CUSTOM_REAL) :: A0, A1, A2, A3, A4, A5, A6, A7, A8
logical :: PML_BOUNDARY_CONDITIONS
ifirstelem = 1
ilastelem = nspec
+ if( PML_BOUNDARY_CONDITIONS ) then
+ potential_dot_dot_acoustic_PML = 0._CUSTOM_REAL
+ PML_dux_dxl = 0._CUSTOM_REAL
+ PML_dux_dzl = 0._CUSTOM_REAL
+ PML_dux_dxl_new = 0._CUSTOM_REAL
+ PML_dux_dzl_new = 0._CUSTOM_REAL
+ endif
+
! loop over spectral elements
do ispec = ifirstelem,ilastelem
@@ -188,7 +196,6 @@
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
@@ -222,38 +229,44 @@
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
+ .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+!------------------------------------------------------------------------------
+!---------------------------- LEFT & RIGHT ------------------------------------
+!------------------------------------------------------------------------------
+ !---------------------- A8 and A6 --------------------------
+ A8 = - d_x_store(iPML) / (k_x_store(iPML) ** 2)
+ A6 = 0.d0
- 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
+ bb = d_x_store(iPML) / k_x_store(iPML) + alpha_x_store(iPML)
+ coef0_x = exp(-bb * deltat)
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = (1.d0 - exp(-bb * deltat / 2.d0)) * A8 / bb
+ coef2_x = (1.d0 - exp(-bb* deltat / 2.d0)) * exp(-bb * deltat / 2.d0) * A8 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A8
+ coef2_x = deltat * A8
+ end if
- 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
+ !---------------------- A5 and A7 --------------------------
+ A5 = d_x_store(iPML)
+ A7 = 0.d0
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+
+ if ( abs( bb ) > 1.d-3) then
+ coef1_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * A5 / bb
+ coef2_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) * &
+ A5 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A5
+ coef2_x = deltat * A5 ! Rene Matzen
+ end if
+
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
@@ -265,147 +278,150 @@
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(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
+ 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
- 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
+ !---------------------------- 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)
- rmemory_acoustic_dux_dx(2,i,j,ispec_PML) = coef0_z*rmemory_acoustic_dux_dx(2,i,j,ispec_PML) &
- + PML_dux_dxl_new(i,j,ispec_PML) * coef1_z + PML_dux_dxl(i,j,ispec_PML) * coef2_z
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_z = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * A8 / bb
+ coef2_z = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * exp(- bb * deltat / 2.d0) * A8 / bb
+ else
+ coef1_z = deltat / 2.0d0 * A8
+ coef2_z = deltat * A8 ! Rene Matzen
+ end if
-
- 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(2,i,j,ispec_PML) = coef0_z*rmemory_acoustic_dux_dx(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(2,i,j,ispec_PML) = coef0_z*rmemory_acoustic_dux_dx(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
+ !---------------------------- 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(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
+ 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
- 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
+ !---------------------------- A6 ----------------------------
+ if ( abs( alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A6 = 0.d0
+ elseif ( abs( alpha_x_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) ) < 1.d-3 &
+ ) then
+ A6 = - d_z_store(iPML) / ( k_z_store(iPML) ** 2 )
+ elseif ( abs( alpha_x_store(iPML) - alpha_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( d_x_store(iPML) - d_z_store(iPML) ) < 1.d-3 .AND. &
+ abs( k_x_store(iPML) - k_z_store(iPML) ) < 1.d-3 &
+ ) then
+ A6 = 0.d0
+ else
+ A6 = d_z_store(iPML) * &
+ ( k_z_store(iPML) * k_x_store(iPML) * ( alpha_x_store(iPML) - alpha_z_store(iPML) ) &
+ + k_z_store(iPML) * d_x_store(iPML) - k_x_store(iPML) * d_z_store(iPML) &
+ ) &
+ / ( k_z_store(iPML)**2 * ( k_z_store(iPML) * ( alpha_z_store(iPML) - alpha_x_store(iPML) ) &
+ + d_z_store(iPML) ) &
+ )
+ end if
+ bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
+ coef0_z = exp(- bb * deltat)
- rmemory_acoustic_dux_dz(2,i,j,ispec_PML) = coef0_z * rmemory_acoustic_dux_dz(2,i,j,ispec_PML) &
- + PML_dux_dzl_new(i,j,ispec_PML) *coef1_z + PML_dux_dzl(i,j,ispec_PML) * coef2_z
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_z = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * A6 / bb
+ coef2_z = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * exp(- bb * deltat / 2.d0) * A6 / bb
+ else
+ coef1_z = deltat / 2.0d0 * A6
+ coef2_z = deltat * A6 ! Rene Matzen
+ end if
- 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(2,i,j,ispec_PML) = coef0_z * rmemory_acoustic_dux_dz(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(2,i,j,ispec_PML) = coef0_z * rmemory_acoustic_dux_dz(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(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) + &
@@ -413,35 +429,42 @@
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
+!------------------------------------------------------------------------------
+!---------------------------- TOP & BOTTOM ------------------------------------
+!------------------------------------------------------------------------------
+ !---------------------- A5 and A7 --------------------------
+ A7 = d_z_store(iPML)
+ A5 = 0.d0
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ if ( abs( bb ) > 1.d-3) then
+ coef1_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * A7 / bb
+ coef2_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) * &
+ A7 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A7
+ coef2_x = deltat * A7 ! Rene Matzen
+ end if
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
+ !---------------------- A8 and A6 --------------------------
+ A6 = - d_z_store(iPML) / ( k_z_store(iPML) ** 2 )
+ A8 = 0.d0
+ bb = d_z_store(iPML) / k_z_store(iPML) + alpha_z_store(iPML)
+ coef0_x = exp(-bb * deltat)
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = (1.d0 - exp(-bb * deltat / 2.d0)) * A6 / bb
+ coef2_x = (1.d0 - exp(-bb* deltat / 2.d0)) * exp(-bb * deltat / 2.d0) * A6 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A6
+ coef2_x = deltat * A6
+ end if
+
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
@@ -485,13 +508,25 @@
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 ------------------------------------
+!------------------------------------------------------------------------------
- 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))
+ !---------------------- A3 and A4 --------------------------
+ A3 = d_x_store(iPML ) * alpha_x_store(iPML) ** 2
+ A4 = 0.d0
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ if ( abs( bb ) > 1.d-3) then
+ coef1_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * A3 / bb
+ coef2_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) * &
+ A3 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A3
+ coef2_x = deltat * A3 ! Rene Matzen
+ end if
+
rmemory_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
@@ -501,67 +536,99 @@
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
- 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
+!------------------------------------------------------------------------------
+!-------------------------------- CORNER --------------------------------------
+!------------------------------------------------------------------------------
- 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
- 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
+ !---------------------------- 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
+ bb = alpha_x_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = ( 1 - exp(- bb * deltat / 2) ) * A3 / bb
+ coef2_x = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) * A3 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A3
+ coef2_x = deltat * A3 ! Rene Matzen
+ end if
+
rmemory_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)+deltat*potential_dot_acoustic(iglob)) * coef1_x &
+ potential_acoustic(iglob) * coef2_x
- endif
+ !---------------------------- 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
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
- 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
+ if ( abs(bb) > 1.d-3 ) then
+ coef1_x = ( 1 - exp(- bb * deltat / 2) ) * A4 / bb
+ coef2_x = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) * A4 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A4
+ coef2_x = deltat * A4 ! Rene Matzen
+ end if
- rmemory_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
-
- 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(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
+
else
- 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 )
+!------------------------------------------------------------------------------
+!-------------------------------- TOP & BOTTOM --------------------------------
+!------------------------------------------------------------------------------
+ !---------------------- A3 and A4 ----------------------------
+ A4 = d_z_store(iPML ) * alpha_z_store(iPML) ** 2
+ A3 = 0.d0
+ bb = alpha_z_store(iPML)
+ coef0_x = exp(- bb * deltat)
+ if ( abs( bb ) > 1.d-3) then
+ coef1_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * A4 / bb
+ coef2_x = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) * &
+ A4 / bb
+ else
+ coef1_x = deltat / 2.0d0 * A4
+ coef2_x = deltat * A4 ! Rene Matzen
+ end if
+
+
rmemory_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 &
@@ -572,30 +639,50 @@
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 ------------------------------------
+!------------------------------------------------------------------------------
+ A0 = - alpha_x_store( iPML ) * d_x_store(iPML)
+ A1 = d_x_store(iPML)
+ A2 = k_x_store(iPML)
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
- (d_x_store(iPML)*potential_dot_acoustic(iglob)+ &
+ (A1*potential_dot_acoustic(iglob)+ &
rmemory_potential_acoustic(1,i,j,ispec_PML)+rmemory_potential_acoustic(2,i,j,ispec_PML)&
- -d_x_store(iPML)*alpha_x_store(iPML)*potential_acoustic(iglob))
+ +A0*potential_acoustic(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
+!------------------------------------------------------------------------------
+!-------------------------------- 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)
+ A1 = d_x_store(iPML) * k_z_store(iPML) + d_z_store(iPML) * k_x_store(iPML)
+
+ A2 = k_x_store(iPML) * k_z_store(iPML)
+
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)+d_z_store(iPML)*k_x_store(iPML))*potential_dot_acoustic(iglob)+ &
+ (A1*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)-&
- 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))
+ +A0*potential_acoustic(iglob))
else
+!------------------------------------------------------------------------------
+!-------------------------------- TOP & BOTTOM --------------------------------
+!------------------------------------------------------------------------------
+ A0 = - alpha_z_store( iPML ) * d_z_store(iPML)
+ A1 = d_z_store(iPML)
+ A2 = k_z_store(iPML)
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
- (d_z_store(iPML)*potential_dot_acoustic(iglob)+ &
+ (A1*potential_dot_acoustic(iglob)+ &
rmemory_potential_acoustic(1,i,j,ispec_PML)+rmemory_potential_acoustic(2,i,j,ispec_PML)&
- -d_z_store(iPML)*alpha_z_store(iPML)*potential_acoustic(iglob))
+ +A0*potential_acoustic(iglob))
endif
endif
@@ -623,18 +710,14 @@
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. &
+ 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) &
- potential_dot_dot_acoustic_PML(i,j,ispec_PML)
- elseif ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+ 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
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(i,j,ispec_PML)
-
else
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- potential_dot_dot_acoustic_PML(i,j,ispec_PML)
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-07-11 21:26:25 UTC (rev 20521)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-07-13 08:58:41 UTC (rev 20522)
@@ -2859,13 +2859,17 @@
!elastic PML memory variables
if (any_elastic .and. nspec_PML>0) then
- allocate(rmemory_displ_elastic(2,3,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_displ_elastic(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)
+ 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)
+ 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)
+ 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)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_duz_dz'
- allocate(rmemory_dux_dx(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_dux_dz(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_duz_dx(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_duz_dz(2,NGLLX,NGLLZ,nspec_PML))
-
rmemory_displ_elastic(:,:,:,:,:) = ZERO
rmemory_dux_dx(:,:,:,:) = ZERO
rmemory_dux_dz(:,:,:,:) = ZERO
@@ -2875,7 +2879,6 @@
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))
@@ -2885,9 +2888,12 @@
if (any_acoustic .and. nspec_PML>0) then
- allocate(rmemory_potential_acoustic(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_acoustic_dux_dx(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_acoustic_dux_dz(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_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)
+ 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)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_acoustic_dux_dz'
rmemory_potential_acoustic = ZERO
rmemory_acoustic_dux_dx = ZERO
@@ -2896,7 +2902,6 @@
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))
@@ -2911,11 +2916,8 @@
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_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))
More information about the CIG-COMMITS
mailing list