[cig-commits] r21934 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Thu Apr 25 05:25:26 PDT 2013
Author: xie.zhinan
Date: 2013-04-25 05:25:25 -0700 (Thu, 25 Apr 2013)
New Revision: 21934
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
Log:
clean the code of CPML
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-04-25 10:17:05 UTC (rev 21933)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-04-25 12:25:25 UTC (rev 21934)
@@ -732,11 +732,11 @@
! because array is_CPML() is unallocated when PML_CONDITIONS is false
if(is_CPML(ispec)) then
! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
- call pml_compute_memory_variables_elastic(ispec,ispec_CPML,deltat,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
- tempx3,tempy3,tempz3,NSPEC_AB,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian)
+ call pml_compute_memory_variables_elastic(ispec,ispec_CPML,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
+ tempx3,tempy3,tempz3)
! calculates contribution from each C-PML element to update acceleration
- call pml_compute_accel_contribution_elastic(ispec,ispec_CPML,deltat,nspec_AB,jacobian)
+ call pml_compute_accel_contribution_elastic(ispec,ispec_CPML)
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-04-25 10:17:05 UTC (rev 21933)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-04-25 12:25:25 UTC (rev 21934)
@@ -26,7 +26,7 @@
!
! United States and French Government Sponsorship Acknowledged.
-subroutine pml_compute_accel_contribution_elastic(ispec,ispec_CPML,deltat,nspec_AB,jacobian)
+subroutine pml_compute_accel_contribution_elastic(ispec,ispec_CPML)
! calculates contribution from each C-PML element to update acceleration to the global mesh
@@ -35,8 +35,8 @@
! Anisotropic-Medium PML for Vector FETD With Modified Basis Functions,
! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006)
- use specfem_par, only: ibool,wgllwgll_yz,wgllwgll_xz,wgllwgll_xy,it,kappastore,rhostore
- use specfem_par_elastic, only: displ,veloc,ispec_is_elastic
+ use specfem_par, only: ibool,wgllwgll_yz,wgllwgll_xz,wgllwgll_xy,jacobian,it,deltat,kappastore,rhostore
+ use specfem_par_elastic, only: displ,veloc
use pml_par, only: NSPEC_CPML,CPML_regions,spec_to_CPML, &
d_store_x,d_store_y,d_store_z,K_store_x,K_store_y,K_store_z,alpha_store,&
rmemory_displ_elastic,accel_elastic_CPML
@@ -45,19 +45,13 @@
implicit none
- integer, intent(in) :: ispec,ispec_CPML,nspec_AB
+ integer, intent(in) :: ispec,ispec_CPML
- real(kind=CUSTOM_REAL), intent(in) :: deltat
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_AB), intent(in) :: jacobian
-
! local parameters
integer :: i,j,k,iglob
-
real(kind=CUSTOM_REAL) :: fac1,fac2,fac3,fac4,rhol,jacobianl
real(kind=CUSTOM_REAL) :: bb,coef0_1,coef1_1,coef2_1,coef0_2,coef1_2,coef2_2,coef0_3,coef1_3,coef2_3
- real(kind=CUSTOM_REAL) :: A0,A1,A2,A3,A4,A5 ! for convolution of acceleration
- real(kind=CUSTOM_REAL) :: temp_A3,temp_A4,temp_A5
+ real(kind=CUSTOM_REAL) :: A0,A1,A2,A3,A4,A5,temp_A3! for convolution of acceleration
do k=1,NGLLZ
do j=1,NGLLY
@@ -65,6 +59,10 @@
rhol = rhostore(i,j,k,ispec)
jacobianl = jacobian(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+ fac4 = sqrt(fac1 * fac2 * fac3)
if( CPML_regions(ispec_CPML) == CPML_X_ONLY ) then
!------------------------------------------------------------------------------
@@ -72,7 +70,6 @@
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -106,39 +103,12 @@
A4 = 0.d0
A5 = 0.d0
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
- fac4 = sqrt(fac1 * fac2 * fac3)
-
- accel_elastic_CPML(1,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(1,iglob) + A2 * displ(1,iglob) + &
- A3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(2,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(2,iglob) + A2 * displ(2,iglob) + &
- A3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(3,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(3,iglob) + A2 * displ(3,iglob) + &
- A3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- Y-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -172,39 +142,12 @@
A4 = 0.d0
A5 = 0.d0
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
- fac4 = sqrt(fac1 * fac2 * fac3)
-
- accel_elastic_CPML(1,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(1,iglob) + A2 * displ(1,iglob) + &
- A3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(2,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(2,iglob) + A2 * displ(2,iglob) + &
- A3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(3,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(3,iglob) + A2 * displ(3,iglob) + &
- A3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- Z-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -238,39 +181,12 @@
A4 = 0.d0
A5 = 0.d0
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
- fac4 = sqrt(fac1 * fac2 * fac3)
-
- accel_elastic_CPML(1,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(1,iglob) + A2 * displ(1,iglob) + &
- A3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(2,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(2,iglob) + A2 * displ(2,iglob) + &
- A3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(3,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(3,iglob) + A2 * displ(3,iglob) + &
- A3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- XY-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -288,22 +204,22 @@
rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) &
+ (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + displ(1,iglob) * coef2_1
rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) &
- + (displ(1,iglob) + deltat * veloc(1,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(1,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * it*deltat * coef1_2 &
+ + displ(1,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) = 0.0
rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) &
+ (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + displ(2,iglob) * coef2_1
rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) &
- + (displ(2,iglob) + deltat * veloc(2,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(2,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * it*deltat * coef1_2 &
+ + displ(2,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) = 0.0
rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) &
+ (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + displ(3,iglob) * coef2_1
rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) &
- + (displ(3,iglob) + deltat * veloc(3,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(3,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * it*deltat * coef1_2 &
+ + displ(3,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) = 0.0
!---------------------- A0, A1, A2, A3, A4 and A5 --------------------------
@@ -316,43 +232,16 @@
A3 = -2.0 * alpha_store(i,j,k,ispec_CPML) * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) &
+ alpha_store(i,j,k,ispec_CPML)**2 * ( d_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
+ d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) ) + alpha_store(i,j,k,ispec_CPML)**2 &
- * (it+0.0) * deltat * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+ * it*deltat * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
A4 = -alpha_store(i,j,k,ispec_CPML)**2 * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
A5 = 0.0
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
- fac4 = sqrt(fac1 * fac2 * fac3)
-
- accel_elastic_CPML(1,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(1,iglob) + A2 * displ(1,iglob) + &
- A3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(2,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(2,iglob) + A2 * displ(2,iglob) + &
- A3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(3,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(3,iglob) + A2 * displ(3,iglob) + &
- A3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- XZ-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -370,22 +259,22 @@
rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) &
+ (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + displ(1,iglob) * coef2_1
rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) &
- + (displ(1,iglob) + deltat * veloc(1,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(1,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * it*deltat * coef1_2 &
+ + displ(1,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(1,i,j,k,ispec_CPML,3)= 0.0
rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) &
+ (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + displ(2,iglob) * coef2_1
rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) &
- + (displ(2,iglob) + deltat * veloc(2,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(2,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * it*deltat * coef1_2 &
+ + displ(2,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(2,i,j,k,ispec_CPML,3)= 0.0
rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) &
+ (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + displ(3,iglob) * coef2_1
rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) &
- + (displ(3,iglob) + deltat * veloc(3,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(3,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * it*deltat * coef1_2 &
+ + displ(3,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(3,i,j,k,ispec_CPML,3)= 0.0
!---------------------- A0, A1, A2, A3, A4 and A5 --------------------------
@@ -395,48 +284,19 @@
A2 = d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) &
- alpha_store(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
- alpha_store(i,j,k,ispec_CPML) * d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML)
- if( ispec_is_elastic(ispec) ) then
- A3 = -2.0 * alpha_store(i,j,k,ispec_CPML) * d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) &
- + alpha_store(i,j,k,ispec_CPML)**2 * ( d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
- + d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) ) + alpha_store(i,j,k,ispec_CPML)**2 &
- * (it+0.0) * deltat * d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
- endif
+ A3 = -2.0 * alpha_store(i,j,k,ispec_CPML) * d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) &
+ + alpha_store(i,j,k,ispec_CPML)**2 * ( d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
+ + d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) ) + alpha_store(i,j,k,ispec_CPML)**2 &
+ * it*deltat * d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
A4 = -alpha_store(i,j,k,ispec_CPML)**2 * d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
A5 = 0.0
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
- fac4 = sqrt(fac1 * fac2 * fac3)
-
- accel_elastic_CPML(1,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(1,iglob) + A2 * displ(1,iglob) + &
- A3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(2,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(2,iglob) + A2 * displ(2,iglob) + &
- A3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(3,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(3,iglob) + A2 * displ(3,iglob) + &
- A3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- YZ-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -454,22 +314,22 @@
rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) &
+ (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + displ(1,iglob) * coef2_1
rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) &
- + (displ(1,iglob) + deltat * veloc(1,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(1,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * it*deltat * coef1_2 &
+ + displ(1,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(1,i,j,k,ispec_CPML,3)=0.d0
rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) &
+ (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + displ(2,iglob) * coef2_1
rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) &
- + (displ(2,iglob) + deltat * veloc(2,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(2,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * it*deltat * coef1_2 &
+ + displ(2,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(2,i,j,k,ispec_CPML,3)=0.d0
rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) &
+ (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + displ(3,iglob) * coef2_1
rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) &
- + (displ(3,iglob) + deltat * veloc(3,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(3,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * it*deltat * coef1_2 &
+ + displ(3,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(3,i,j,k,ispec_CPML,3)=0.d0
!---------------------- A0, A1, A2, A3, A4 and A5 --------------------------
@@ -482,43 +342,16 @@
A3 = -2.0 * alpha_store(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) &
+ alpha_store(i,j,k,ispec_CPML)**2 * ( d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
+ d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) ) + alpha_store(i,j,k,ispec_CPML)**2 &
- * (it+0.0) * deltat * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+ * it*deltat * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
A4 = -alpha_store(i,j,k,ispec_CPML)**2 * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
A5 = 0.0
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
- fac4 = sqrt(fac1 * fac2 * fac3)
-
- accel_elastic_CPML(1,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(1,iglob) + A2 * displ(1,iglob) + &
- A3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(2,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(2,iglob) + A2 * displ(2,iglob) + &
- A3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
- )
-
- accel_elastic_CPML(3,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(3,iglob) + A2 * displ(3,iglob) + &
- A3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
!------------------------------------------------------------------------------
!---------------------------- XYZ-corner C-PML --------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -540,29 +373,29 @@
rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) &
+ (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + displ(1,iglob) * coef2_1
rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) &
- + (displ(1,iglob) + deltat * veloc(1,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(1,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * it*deltat * coef1_2 &
+ + displ(1,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) = coef0_3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
- + (displ(1,iglob) + deltat * veloc(1,iglob)) * ((it+0.0) * deltat)**2 * coef1_3 &
- + displ(1,iglob) * ((it-0.0) * deltat)**2 * coef2_3
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * (it*deltat)**2 * coef1_3 &
+ + displ(1,iglob) * (it*deltat)**2 * coef2_3
rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) &
+ (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + displ(2,iglob) * coef2_1
rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) &
- + (displ(2,iglob) + deltat * veloc(2,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(2,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * it*deltat * coef1_2 &
+ + displ(2,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) = coef0_3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
- + (displ(2,iglob) + deltat * veloc(2,iglob)) * ((it+0.0) * deltat)**2 * coef1_3 &
- + displ(2,iglob) * ((it-0.0) * deltat)**2 * coef2_3
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * (it*deltat)**2 * coef1_3 &
+ + displ(2,iglob) * (it*deltat)**2 * coef2_3
rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) = coef0_1 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) &
+ (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + displ(3,iglob) * coef2_1
rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) = coef0_2 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) &
- + (displ(3,iglob) + deltat * veloc(3,iglob)) * (it+0.0) * deltat * coef1_2 &
- + displ(3,iglob) * (it-0.0) * deltat * coef2_2
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * it*deltat * coef1_2 &
+ + displ(3,iglob) * it*deltat * coef2_2
rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) = coef0_3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
- + (displ(3,iglob) + deltat * veloc(3,iglob)) * ((it+0.0) * deltat)**2 * coef1_3 &
- + displ(3,iglob) * ((it-0.0) * deltat)**2 * coef2_3
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * (it*deltat)**2 * coef1_3 &
+ + displ(3,iglob) * (it*deltat)**2 * coef2_3
!---------------------- A0, A1, A2, A3, A4 and A5 --------------------------
A0 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML)
@@ -586,18 +419,16 @@
d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) + &
d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
)
- temp_A4 = -2.0 * alpha_store(i,j,k,ispec_CPML) * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * &
- d_store_z(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)**2 * ( &
- d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) + &
- d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) + &
- d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
- )
- temp_A5 = 0.5 * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
-
+! temp_A4 = -2.0 * alpha_store(i,j,k,ispec_CPML) * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * &
+! d_store_z(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)**2 * ( &
+! d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) + &
+! d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) + &
+! d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+! )
+! temp_A5 = 0.5 * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
! A3 = temp_A3 + (it+0.0)*deltat*temp_A4 + ((it+0.0)*deltat)**2*temp_A5
! A4 = -temp_A4 -2.0*(it+0.0)*deltat*temp_A5
! A5 = temp_A5
-
!ZN the full experssion of A3,A4,A5 are given by above equation, here we use reduced
!ZN exprssion of A3,A4,A5 in order to stabilized the code.
@@ -605,33 +436,28 @@
A4 = 0.0
A5 = 0.0
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
- fac4 = sqrt(fac1 * fac2 * fac3)
+ endif
- accel_elastic_CPML(1,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(1,iglob) + A2 * displ(1,iglob) + &
- A3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
- )
+ accel_elastic_CPML(1,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
+ ( A1 * veloc(1,iglob) + A2 * displ(1,iglob) + &
+ A3 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,1) + &
+ A4 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,2) + &
+ A5 * rmemory_displ_elastic(1,i,j,k,ispec_CPML,3) &
+ )
- accel_elastic_CPML(2,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(2,iglob) + A2 * displ(2,iglob) + &
- A3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
- )
+ accel_elastic_CPML(2,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
+ ( A1 * veloc(2,iglob) + A2 * displ(2,iglob) + &
+ A3 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,1) + &
+ A4 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,2) + &
+ A5 * rmemory_displ_elastic(2,i,j,k,ispec_CPML,3) &
+ )
- accel_elastic_CPML(3,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
- ( A1 * veloc(3,iglob) + A2 * displ(3,iglob) + &
- A3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) + &
- A4 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) + &
- A5 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
- )
-
- endif
+ accel_elastic_CPML(3,i,j,k,ispec_CPML) = fac4 * rhol * jacobianl * &
+ ( A1 * veloc(3,iglob) + A2 * displ(3,iglob) + &
+ A3 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,1) + &
+ A4 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,2) + &
+ A5 * rmemory_displ_elastic(3,i,j,k,ispec_CPML,3) &
+ )
enddo
enddo
enddo
@@ -668,8 +494,7 @@
integer :: i,j,k,iglob
real(kind=CUSTOM_REAL) :: fac1,fac2,fac3,fac4,kappal,jacobianl
real(kind=CUSTOM_REAL) :: bb,coef0_1,coef1_1,coef2_1,coef0_2,coef1_2,coef2_2,coef0_3,coef1_3,coef2_3
- real(kind=CUSTOM_REAL) :: A0,A1,A2,A3,A4,A5 ! for convolution of acceleration
- real(kind=CUSTOM_REAL) :: temp_A3,temp_A4,temp_A5
+ real(kind=CUSTOM_REAL) :: A0,A1,A2,A3,A4,A5,temp_A3 ! for convolution of acceleration
do k=1,NGLLZ
do j=1,NGLLY
@@ -688,7 +513,6 @@
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -713,20 +537,12 @@
A4 = 0.d0
A5 = 0.d0
- potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML) = fac4 * 1.0/kappal *jacobianl * &
- ( A1 * potential_dot_acoustic(iglob) + A2 * potential_acoustic(iglob) + &
- A3 * rmemory_potential_acoustic(i,j,k,ispec_CPML,1)+ &
- A4 * rmemory_potential_acoustic(i,j,k,ispec_CPML,2)+ &
- A5 * rmemory_potential_acoustic(i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- Y-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -751,20 +567,12 @@
A4 = 0.d0
A5 = 0.d0
- potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML) = fac4 * 1.0/kappal *jacobianl * &
- ( A1 * potential_dot_acoustic(iglob) + A2 * potential_acoustic(iglob) + &
- A3 * rmemory_potential_acoustic(i,j,k,ispec_CPML,1)+ &
- A4 * rmemory_potential_acoustic(i,j,k,ispec_CPML,2)+ &
- A5 * rmemory_potential_acoustic(i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- Z-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -789,20 +597,12 @@
A4 = 0.d0
A5 = 0.d0
- potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML) = fac4 * 1.0/kappal *jacobianl * &
- ( A1 * potential_dot_acoustic(iglob) + A2 * potential_acoustic(iglob) + &
- A3 * rmemory_potential_acoustic(i,j,k,ispec_CPML,1)+ &
- A4 * rmemory_potential_acoustic(i,j,k,ispec_CPML,2)+ &
- A5 * rmemory_potential_acoustic(i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- XY-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -839,20 +639,12 @@
A4 = -alpha_store(i,j,k,ispec_CPML)**2 * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
A5 = 0.0
- potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML) = fac4 * 1.0/kappal *jacobianl * &
- ( A1 * potential_dot_acoustic(iglob) + A2 * potential_acoustic(iglob) + &
- A3 * rmemory_potential_acoustic(i,j,k,ispec_CPML,1)+ &
- A4 * rmemory_potential_acoustic(i,j,k,ispec_CPML,2)+ &
- A5 * rmemory_potential_acoustic(i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- XZ-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -889,20 +681,12 @@
A4 = -alpha_store(i,j,k,ispec_CPML)**2 * d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
A5 = 0.0
- potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML) = fac4 * 1.0/kappal *jacobianl * &
- ( A1 * potential_dot_acoustic(iglob) + A2 * potential_acoustic(iglob) + &
- A3 * rmemory_potential_acoustic(i,j,k,ispec_CPML,1)+ &
- A4 * rmemory_potential_acoustic(i,j,k,ispec_CPML,2)+ &
- A5 * rmemory_potential_acoustic(i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- YZ-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -939,20 +723,12 @@
A4 = -alpha_store(i,j,k,ispec_CPML)**2 * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
A5 = 0.0
- potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML) = fac4 * 1.0/kappal *jacobianl * &
- ( A1 * potential_dot_acoustic(iglob) + A2 * potential_acoustic(iglob) + &
- A3 * rmemory_potential_acoustic(i,j,k,ispec_CPML,1)+ &
- A4 * rmemory_potential_acoustic(i,j,k,ispec_CPML,2)+ &
- A5 * rmemory_potential_acoustic(i,j,k,ispec_CPML,3) &
- )
-
else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
!------------------------------------------------------------------------------
!---------------------------- XYZ-corner C-PML --------------------------------
!------------------------------------------------------------------------------
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1003,14 +779,13 @@
d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) + &
d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
)
- temp_A4 = -2.0 * alpha_store(i,j,k,ispec_CPML) * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * &
- d_store_z(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)**2 * ( &
- d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) + &
- d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) + &
- d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
- )
- temp_A5 = 0.5 * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
-
+! temp_A4 = -2.0 * alpha_store(i,j,k,ispec_CPML) * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * &
+! d_store_z(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)**2 * ( &
+! d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) + &
+! d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) + &
+! d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+! )
+! temp_A5 = 0.5 * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
! A3 = temp_A3 + it*deltat*temp_A4 + (it*deltat)**2*temp_A5
! A4 = -temp_A4 -2.0*it*deltat*temp_A5
! A5 = temp_A5
@@ -1022,13 +797,14 @@
A4 = 0.0
A5 = 0.0
- potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML) = fac4 * 1.0/kappal *jacobianl * &
- ( A1 * potential_dot_acoustic(iglob) + A2 * potential_acoustic(iglob) + &
- A3 * rmemory_potential_acoustic(i,j,k,ispec_CPML,1)+ &
- A4 * rmemory_potential_acoustic(i,j,k,ispec_CPML,2)+ &
- A5 * rmemory_potential_acoustic(i,j,k,ispec_CPML,3) &
- )
endif
+
+ potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML) = fac4 * 1.0/kappal *jacobianl * &
+ ( A1 * potential_dot_acoustic(iglob) + A2 * potential_acoustic(iglob) + &
+ A3 * rmemory_potential_acoustic(i,j,k,ispec_CPML,1)+ &
+ A4 * rmemory_potential_acoustic(i,j,k,ispec_CPML,2)+ &
+ A5 * rmemory_potential_acoustic(i,j,k,ispec_CPML,3) &
+ )
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-04-25 10:17:05 UTC (rev 21933)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-04-25 12:25:25 UTC (rev 21934)
@@ -26,9 +26,8 @@
!
! United States and French Government Sponsorship Acknowledged.
-subroutine pml_compute_memory_variables_elastic(ispec,ispec_CPML,deltat,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
- tempx3,tempy3,tempz3,NSPEC_AB,xix,xiy,xiz, &
- etax,etay,etaz,gammax,gammay,gammaz,jacobian)
+subroutine pml_compute_memory_variables_elastic(ispec,ispec_CPML,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
+ tempx3,tempy3,tempz3)
! calculates C-PML elastic memory variables and computes stress sigma
! second-order accurate convolution term calculation from equation (21) of
@@ -36,28 +35,23 @@
! Anisotropic-Medium PML for Vector FETD With Modified Basis Functions,
! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006)
- use specfem_par, only: it,kappastore,mustore,rhostore,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz
+ use specfem_par, only: wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,it,deltat, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian, &
+ kappastore,mustore,rhostore
use pml_par
use constants, only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS, &
CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ
implicit none
- integer, intent(in) :: ispec,ispec_CPML,NSPEC_AB
-
- real(kind=CUSTOM_REAL), intent(in) :: deltat
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB), intent(in) :: &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian
- real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
-
-
+ integer, intent(in) :: ispec,ispec_CPML
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ), intent(out) :: tempx1,tempx2,tempx3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ), intent(out) :: tempy1,tempy2,tempy3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ), intent(out) :: tempz1,tempz2,tempz3
! local parameters
integer :: i,j,k
-
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
real(kind=CUSTOM_REAL) :: sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz,sigma_yx,sigma_zx,sigma_zy
real(kind=CUSTOM_REAL) :: lambdal,mul,lambdalplus2mul,kappal
real(kind=CUSTOM_REAL) :: duxdxl_x,duxdyl_x,duxdzl_x,duydxl_x,duydyl_x,duzdxl_x,duzdzl_x
@@ -97,7 +91,6 @@
A8 = - d_store_x(i,j,k,ispec_CPML) / (k_store_x(i,j,k,ispec_CPML)**2)
bb = d_store_x(i,j,k,ispec_CPML) / k_store_x(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -134,7 +127,6 @@
A11 = 0.d0
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -171,7 +163,6 @@
A14 = 0.d0
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -201,128 +192,101 @@
duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
+ A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
- !---------------------- A15 and A16 --------------------------
- A15 = k_store_x(i,j,k,ispec_CPML)
- A16 = d_store_x(i,j,k,ispec_CPML)
+ !---------------------- A15 and A16 --------------------------
+ A15 = k_store_x(i,j,k,ispec_CPML)
+ A16 = d_store_x(i,j,k,ispec_CPML)
- bb = alpha_store(i,j,k,ispec_CPML)
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
- coef0_1 = exp(-bb * deltat)
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
- if( abs(bb) > 1.d-5 ) then
- coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
- coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
- else
- coef1_1 = deltat/2.0d0
- coef2_1 = deltat/2.0d0
- endif
+ rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) &
+ + PML_duz_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dzl(i,j,k,ispec_CPML) * coef2_1
+ rmemory_duz_dzl_y(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) &
- + PML_duz_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dzl(i,j,k,ispec_CPML) * coef2_1
- rmemory_duz_dzl_y(i,j,k,ispec_CPML,2) = 0.d0
+ rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) &
+ + PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_1
+ rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) &
- + PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_1
- rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
+ duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
+ + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
+ duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
+ + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
- duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
- duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
+ rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) &
+ + PML_duy_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dzl(i,j,k,ispec_CPML) * coef2_1
+ rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) &
- + PML_duy_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dzl(i,j,k,ispec_CPML) * coef2_1
- rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
+ rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) &
+ + PML_duy_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dyl(i,j,k,ispec_CPML) * coef2_1
+ rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) &
- + PML_duy_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dyl(i,j,k,ispec_CPML) * coef2_1
- rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
+ duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
+ + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
+ duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
+ + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
- duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
- duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
+ !---------------------- A17 and A18 --------------------------
+ A17 = 1.d0
+ A18 = 0.0
- !---------------------- A17 and A18 --------------------------
- A17 = 1.d0
- A18 = 0.0
+ rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) = 0.d0
+ rmemory_duz_dzl_x(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) = 0.d0
- rmemory_duz_dzl_x(i,j,k,ispec_CPML,2) = 0.d0
+ rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
+ rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
- rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
+ duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
+ + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
+ duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
+ + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
- duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
- duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
+ rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = 0.d0
+ rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = 0.d0
- rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
+ rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) = 0.d0
+ rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) = 0.d0
- rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
+ duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
+ + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
+ duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
- duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
- duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
+ !---------------------- A19 and A20 --------------------------
+ A19 = 1.d0
+ A20 = 0.0
- !---------------------- A19 and A20 --------------------------
- A19 = 1.d0
- A20 = 0.0
+ rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) = 0.d0
+ rmemory_duy_dyl_x(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) = 0.d0
- rmemory_duy_dyl_x(i,j,k,ispec_CPML,2) = 0.d0
+ rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
+ rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
- rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
+ duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
+ + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
+ duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
+ + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
- duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
- duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
+ rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = 0.d0
+ rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = 0.d0
- rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
+ rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) = 0.d0
+ rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
- rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) = 0.d0
- rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
+ duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
+ + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
+ duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
- duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
- duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
- sigma_yx = mul*duxdyl_x + mul*duydxl_x
- sigma_zx = mul*duzdxl_x + mul*duxdzl_x
-
- sigma_xy = mul*duxdyl_y + mul*duydxl_y
- sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
- sigma_zy = mul*duzdyl_y + mul*duydzl_y
-
- sigma_xz = mul*duzdxl_z + mul*duxdzl_z
- sigma_yz = mul*duzdyl_z + mul*duydzl_z
- sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
-
- ! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
-
else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- Y-surface C-PML ---------------------------------
@@ -334,7 +298,6 @@
A8 = 0.d0
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -370,7 +333,6 @@
A11 = - d_store_y(i,j,k,ispec_CPML) / (k_store_y(i,j,k,ispec_CPML) ** 2)
bb = d_store_y(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -406,7 +368,6 @@
A14 = 0.d0
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -467,7 +428,6 @@
A18 = d_store_y(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -531,32 +491,6 @@
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
- sigma_yx = mul*duxdyl_x + mul*duydxl_x
- sigma_zx = mul*duzdxl_x + mul*duxdzl_x
-
- sigma_xy = mul*duxdyl_y + mul*duydxl_y
- sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
- sigma_zy = mul*duzdyl_y + mul*duydzl_y
-
- sigma_xz = mul*duzdxl_z + mul*duxdzl_z
- sigma_yz = mul*duzdyl_z + mul*duydzl_z
- sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
-
- ! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
-
else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
!------------------------------------------------------------------------------
@@ -569,7 +503,6 @@
A8 = 0.d0
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -605,7 +538,6 @@
A11 = 0.d0
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -641,7 +573,6 @@
A14 = - d_store_z(i,j,k,ispec_CPML) / (k_store_z(i,j,k,ispec_CPML) ** 2)
bb = d_store_z(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -728,7 +659,6 @@
A20 = d_store_z(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -765,32 +695,6 @@
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
- sigma_yx = mul*duxdyl_x + mul*duydxl_x
- sigma_zx = mul*duzdxl_x + mul*duxdzl_x
-
- sigma_xy = mul*duxdyl_y + mul*duydxl_y
- sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
- sigma_zy = mul*duzdyl_y + mul*duydzl_y
-
- sigma_xz = mul*duzdxl_z + mul*duxdzl_z
- sigma_yz = mul*duzdyl_z + mul*duydzl_z
- sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
-
- ! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
-
else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
!------------------------------------------------------------------------------
@@ -804,7 +708,6 @@
d_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) ) / k_store_x(i,j,k,ispec_CPML)**2
bb = d_store_x(i,j,k,ispec_CPML) / k_store_x(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -842,7 +745,6 @@
d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) ) / k_store_y(i,j,k,ispec_CPML)**2
bb = d_store_y(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -880,7 +782,6 @@
A14 = - d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -925,7 +826,6 @@
A16 = d_store_x(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -967,7 +867,6 @@
A18 = d_store_y(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1030,32 +929,6 @@
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
- sigma_yx = mul*duxdyl_x + mul*duydxl_x
- sigma_zx = mul*duzdxl_x + mul*duxdzl_x
-
- sigma_xy = mul*duxdyl_y + mul*duydxl_y
- sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
- sigma_zy = mul*duzdyl_y + mul*duydzl_y
-
- sigma_xz = mul*duzdxl_z + mul*duxdzl_z
- sigma_yz = mul*duzdyl_z + mul*duydzl_z
- sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
-
- ! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
-
else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
!------------------------------------------------------------------------------
@@ -1069,7 +942,6 @@
d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) ) / k_store_x(i,j,k,ispec_CPML)**2
bb = d_store_x(i,j,k,ispec_CPML) / k_store_x(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1107,7 +979,6 @@
A11 = - d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1154,7 +1025,6 @@
- d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) ) / k_store_z(i,j,k,ispec_CPML)**2
bb = d_store_z(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1189,7 +1059,6 @@
A16 = d_store_x(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1257,7 +1126,6 @@
A20 = d_store_z(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1294,32 +1162,6 @@
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
- sigma_yx = mul*duxdyl_x + mul*duydxl_x
- sigma_zx = mul*duzdxl_x + mul*duxdzl_x
-
- sigma_xy = mul*duxdyl_y + mul*duydxl_y
- sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
- sigma_zy = mul*duzdyl_y + mul*duydzl_y
-
- sigma_xz = mul*duzdxl_z + mul*duxdzl_z
- sigma_yz = mul*duzdyl_z + mul*duydzl_z
- sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
-
- ! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
-
else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
!------------------------------------------------------------------------------
@@ -1334,7 +1176,6 @@
A8 = - d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1381,7 +1222,6 @@
d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) ) / k_store_y(i,j,k,ispec_CPML)**2
bb = d_store_y(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1418,7 +1258,6 @@
d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) ) / k_store_z(i,j,k,ispec_CPML)**2
bb = d_store_z(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1479,7 +1318,6 @@
A18 = d_store_y(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1521,7 +1359,6 @@
A20 = d_store_z(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1558,32 +1395,6 @@
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
- sigma_yx = mul*duxdyl_x + mul*duydxl_x
- sigma_zx = mul*duzdxl_x + mul*duxdzl_x
-
- sigma_xy = mul*duxdyl_y + mul*duydxl_y
- sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
- sigma_zy = mul*duzdyl_y + mul*duydzl_y
-
- sigma_xz = mul*duzdxl_z + mul*duxdzl_z
- sigma_yz = mul*duzdyl_z + mul*duydzl_z
- sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
-
- ! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
-
else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
!------------------------------------------------------------------------------
@@ -1608,7 +1419,6 @@
endif
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1620,7 +1430,6 @@
endif
bb = d_store_x(i,j,k,ispec_CPML) / k_store_x(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1683,7 +1492,6 @@
endif
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1695,7 +1503,6 @@
endif
bb = d_store_y(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML) + alpha_store(i,j,k,ispec_CPML)
-
coef0_2 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1757,7 +1564,6 @@
endif
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1818,7 +1624,6 @@
A16 = d_store_x(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1860,7 +1665,6 @@
A18 = d_store_y(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1902,7 +1706,6 @@
A20 = d_store_z(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -1939,34 +1742,36 @@
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
- sigma_yx = mul*duxdyl_x + mul*duydxl_x
- sigma_zx = mul*duzdxl_x + mul*duxdzl_x
+ else
+ stop 'wrong PML flag in PML memory variable calculation routine'
+ endif
- sigma_xy = mul*duxdyl_y + mul*duydxl_y
- sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
- sigma_zy = mul*duzdyl_y + mul*duydzl_y
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
+ sigma_yx = mul*duxdyl_x + mul*duydxl_x
+ sigma_zx = mul*duzdxl_x + mul*duxdzl_x
- sigma_xz = mul*duzdxl_z + mul*duxdzl_z
- sigma_yz = mul*duzdyl_z + mul*duydzl_z
- sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
+ sigma_xy = mul*duxdyl_y + mul*duydxl_y
+ sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
+ sigma_zy = mul*duzdyl_y + mul*duydzl_y
- ! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+ sigma_xz = mul*duzdxl_z + mul*duxdzl_z
+ sigma_yz = mul*duzdyl_z + mul*duydzl_z
+ sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
- else
- stop 'wrong PML flag in PML memory variable calculation routine'
- endif
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+
enddo
enddo
enddo
@@ -2100,10 +1905,6 @@
dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
+ A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
- temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
-
else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- Y-surface C-PML ---------------------------------
@@ -2178,10 +1979,6 @@
dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
+ A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
- temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
-
else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
!------------------------------------------------------------------------------
@@ -2257,10 +2054,6 @@
dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
+ A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
- temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
-
else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
!------------------------------------------------------------------------------
@@ -2346,10 +2139,6 @@
dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
+ A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
- temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
-
else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
!------------------------------------------------------------------------------
@@ -2435,9 +2224,6 @@
dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
+ A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
- temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
@@ -2524,10 +2310,6 @@
dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
+ A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
- temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
-
else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
!------------------------------------------------------------------------------
@@ -2660,7 +2442,6 @@
endif
bb = alpha_store(i,j,k,ispec_CPML)
-
coef0_1 = exp(-bb * deltat)
if( abs(bb) > 1.d-5 ) then
@@ -2697,13 +2478,14 @@
dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
+ A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
- temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
-
else
stop 'wrong PML flag in PML memory variable calculation routine'
endif
+
+ temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
+ temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
+ temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
+
enddo
enddo
enddo
More information about the CIG-COMMITS
mailing list