[cig-commits] r21686 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Fri Mar 29 11:04:32 PDT 2013
Author: xie.zhinan
Date: 2013-03-29 11:04:31 -0700 (Fri, 29 Mar 2013)
New Revision: 21686
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
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:
fix some bugs in code related to PML
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-03-29 15:18:49 UTC (rev 21685)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-03-29 18:04:31 UTC (rev 21686)
@@ -90,7 +90,9 @@
ALPHA_MAX_PML = PI*f0_FOR_PML ! ELASTIC from Festa and Vilotte (2005)
ALPHA_MAX_PML = PI*f0_FOR_PML*2.0 ! ACOUSTIC from Festa and Vilotte (2005)
- CPML_width_x = CPML_width
+ !the idea of fix PML width is bad when element sizes are not the same in x,y,z direction
+
+ CPML_width_x = CPML_width
CPML_width_y = CPML_width
CPML_width_z = CPML_width
@@ -99,12 +101,14 @@
y_min = minval(ystore(:))
y_max = maxval(ystore(:))
z_min = minval(zstore(:))
+ z_max = maxval(zstore(:))
x_min_all = 10.e30
x_max_all = -10.e30
y_min_all = 10.e30
y_max_all = -10.e30
z_min_all = 10.e30
+ z_max_all = -10.e30
call min_all_all_cr(x_min,x_min_all)
call min_all_all_cr(y_min,y_min_all)
@@ -112,28 +116,20 @@
call max_all_all_cr(x_max,x_max_all)
call max_all_all_cr(y_max,y_max_all)
+ call max_all_all_cr(z_max,z_max_all)
+
x_origin = (x_min_all + x_max_all)/2.0
y_origin = (y_min_all + y_max_all)/2.0
- z_origin = z_min_all/2.0
+ z_origin = (z_max_all + z_min_all)/2.0
! determines equations of C-PML/mesh interface planes
-!ZN xoriginleft = minval(xstore(:)) + CPML_width_x
-!ZN xoriginright = maxval(xstore(:)) - CPML_width_x
-!ZN yoriginback = minval(ystore(:)) + CPML_width_y
-!ZN yoriginfront = maxval(ystore(:)) - CPML_width_y
-!ZN zoriginbottom = minval(zstore(:)) + CPML_width_z
-
xoriginleft = x_min_all + CPML_width_x
xoriginright = x_max_all - CPML_width_x
yoriginback = y_min_all + CPML_width_y
yoriginfront = y_max_all - CPML_width_y
zoriginbottom = z_min_all + CPML_width_z
- if( PML_INSTEAD_OF_FREE_SURFACE ) then
-!ZN zorigintop = maxval(zstore(:)) - CPML_width_z
- z_max = maxval(zstore(:))
- z_max_all = -10.e30
- call max_all_all_cr(z_max,z_max_all)
+ if( PML_INSTEAD_OF_FREE_SURFACE ) then
zorigintop = z_max_all - CPML_width_z
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-03-29 15:18:49 UTC (rev 21685)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-03-29 18:04:31 UTC (rev 21686)
@@ -90,7 +90,6 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
- 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) :: hp1,hp2,hp3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: chi_elem
@@ -98,8 +97,6 @@
real(kind=CUSTOM_REAL) :: temp1l,temp2l,temp3l
real(kind=CUSTOM_REAL) :: temp1l_new,temp2l_new,temp3l_new
- real(kind=CUSTOM_REAL) :: lambdal,mul,lambdalplus2mul
-
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
real(kind=CUSTOM_REAL) :: dpotentialdxl,dpotentialdyl,dpotentialdzl
real(kind=CUSTOM_REAL) :: dpotentialdxl_new,dpotentialdyl_new,dpotentialdzl_new
@@ -233,13 +230,11 @@
! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
if(CPML_mask_ibool(ispec)) then
! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
- call pml_compute_memory_variables(ispec,ispec_CPML,deltat,jacobianl,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
- tempx3,tempy3,tempz3,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz, &
- sigma_yx,sigma_zx,sigma_zy,lambdal,mul,lambdalplus2mul,xixl,xiyl,xizl, &
- etaxl,etayl,etazl,gammaxl,gammayl,gammazl)
+ call pml_compute_memory_variables(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)
! calculates contribution from each C-PML element to update acceleration
- call pml_compute_accel_contribution(ispec,ispec_CPML,deltat,jacobianl,accel_elastic_CPML)
+ call pml_compute_accel_contribution(ispec,ispec_CPML,deltat,nspec_AB,jacobian,potential_dot_dot_acoustic_CPML)
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-03-29 15:18:49 UTC (rev 21685)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-03-29 18:04:31 UTC (rev 21686)
@@ -729,13 +729,11 @@
! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
if(CPML_mask_ibool(ispec)) then
! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
- call pml_compute_memory_variables(ispec,ispec_CPML,deltat,jacobianl,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
- tempx3,tempy3,tempz3,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz, &
- sigma_yx,sigma_zx,sigma_zy,lambdal,mul,lambdalplus2mul,NSPEC_AB,xix,xiy,xiz, &
- etax,etay,etaz,gammax,gammay,gammaz,jacobian)
+ call pml_compute_memory_variables(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)
! calculates contribution from each C-PML element to update acceleration
- call pml_compute_accel_contribution(ispec,ispec_CPML,deltat,jacobianl,accel_elastic_CPML)
+ call pml_compute_accel_contribution(ispec,ispec_CPML,deltat,nspec_AB,jacobian,accel_elastic_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-03-29 15:18:49 UTC (rev 21685)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-03-29 18:04:31 UTC (rev 21686)
@@ -26,7 +26,7 @@
!
! United States and French Government Sponsorship Acknowledged.
-subroutine pml_compute_accel_contribution(ispec,ispec_CPML,deltat,jacobianl,accel_elastic_CPML)
+subroutine pml_compute_accel_contribution(ispec,ispec_CPML,deltat,nspec_AB,jacobian,accel_elastic_CPML)
! calculates contribution from each C-PML element to update acceleration to the global mesh
@@ -44,16 +44,17 @@
implicit none
- integer, intent(in) :: ispec,ispec_CPML
+ integer, intent(in) :: ispec,ispec_CPML,nspec_AB
- real(kind=CUSTOM_REAL), intent(in) :: deltat,jacobianl
+ real(kind=CUSTOM_REAL), intent(in) :: deltat
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ,NSPEC_CPML), intent(out) :: accel_elastic_CPML
+ 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,kappal
+ real(kind=CUSTOM_REAL) :: fac1,fac2,fac3,fac4,rhol,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
@@ -63,8 +64,10 @@
do i=1,NGLLX
if( ispec_is_elastic(ispec) ) then
rhol = rhostore(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
else if( ispec_is_acoustic(ispec) ) then
kappal = kappastore(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
else
stop 'CPML elements should be either acoustic or elastic; exiting...'
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-03-29 15:18:49 UTC (rev 21685)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-03-29 18:04:31 UTC (rev 21686)
@@ -26,9 +26,8 @@
!
! United States and French Government Sponsorship Acknowledged.
-subroutine pml_compute_memory_variables(ispec,ispec_CPML,deltat,jacobianl,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
- tempx3,tempy3,tempz3,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz, &
- sigma_yx,sigma_zx,sigma_zy,lambdal,mul,lambdalplus2mul,NSPEC_AB,xix,xiy,xiz, &
+subroutine pml_compute_memory_variables(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)
! calculates C-PML elastic memory variables and computes stress sigma
@@ -37,17 +36,17 @@
! 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
+ use specfem_par, only: it,kappastore,mustore
use specfem_par_elastic, only: ispec_is_elastic
use specfem_par_acoustic, only: ispec_is_acoustic
use pml_par
- use constants, only: NGLLX,NGLLY,NGLLZ
+ use constants, only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS
implicit none
integer, intent(in) :: ispec,ispec_CPML,NSPEC_AB
- real(kind=CUSTOM_REAL), intent(in) :: lambdal,mul,lambdalplus2mul,deltat
+ 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
@@ -61,6 +60,7 @@
integer :: i,j,k
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
real(kind=CUSTOM_REAL) :: duxdxl_y,duxdyl_y,duydxl_y,duydyl_y,duydzl_y,duzdyl_y,duzdzl_y
real(kind=CUSTOM_REAL) :: duxdxl_z,duxdzl_z,duydyl_z,duydzl_z,duzdxl_z,duzdyl_z,duzdzl_z
@@ -69,11 +69,26 @@
real(kind=CUSTOM_REAL) :: A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17 ! for convolution of strain(complex)
real(kind=CUSTOM_REAL) :: A18,A19,A20 ! for convolution of strain(simple)
- if( CPML_regions(ispec_CPML) == 1 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ kappal = kappastore(i,j,k,ispec)
+ mul = mustore(i,j,k,ispec)
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
+ if( CPML_regions(ispec_CPML) == 1 ) then
+
!------------------------------------------------------------------------------
!---------------------------- X-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
@@ -333,16 +348,6 @@
sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
! form dot product with test vector, non-symmetric form
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
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
@@ -355,14 +360,8 @@
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
endif
- enddo
- enddo
- enddo
- else if( CPML_regions(ispec_CPML) == 2 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ else if( CPML_regions(ispec_CPML) == 2 ) then
!------------------------------------------------------------------------------
!---------------------------- Y-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
@@ -621,16 +620,6 @@
sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
! form dot product with test vector, non-symmetric form
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
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
@@ -643,14 +632,9 @@
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
endif
- enddo
- enddo
- enddo
- else if( CPML_regions(ispec_CPML) == 3 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ else if( CPML_regions(ispec_CPML) == 3 ) then
+
!------------------------------------------------------------------------------
!---------------------------- Z-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
@@ -908,16 +892,6 @@
sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
! form dot product with test vector, non-symmetric form
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
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
@@ -930,14 +904,9 @@
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
endif
- enddo
- enddo
- enddo
- else if( CPML_regions(ispec_CPML) == 4 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ else if( CPML_regions(ispec_CPML) == 4 ) then
+
!------------------------------------------------------------------------------
!---------------------------- XY-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
@@ -1234,16 +1203,6 @@
sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
! form dot product with test vector, non-symmetric form
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
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
@@ -1256,14 +1215,8 @@
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
endif
- enddo
- enddo
- enddo
- else if( CPML_regions(ispec_CPML) == 5 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ else if( CPML_regions(ispec_CPML) == 5 ) then
!------------------------------------------------------------------------------
!---------------------------- XZ-edge C-PML -----------------------------------
@@ -1561,16 +1514,6 @@
sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
! form dot product with test vector, non-symmetric form
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
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
@@ -1583,14 +1526,9 @@
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
endif
- enddo
- enddo
- enddo
- else if( CPML_regions(ispec_CPML) == 6 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ else if( CPML_regions(ispec_CPML) == 6 ) then
+
!------------------------------------------------------------------------------
!---------------------------- YZ-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
@@ -1886,16 +1824,6 @@
sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
! form dot product with test vector, non-symmetric form
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
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
@@ -1908,14 +1836,9 @@
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
endif
- enddo
- enddo
- enddo
- else if( CPML_regions(ispec_CPML) == 7 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ else if( CPML_regions(ispec_CPML) == 7 ) then
+
!------------------------------------------------------------------------------
!---------------------------- XYZ-corner C-PML --------------------------------
!------------------------------------------------------------------------------
@@ -2365,16 +2288,6 @@
sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
! form dot product with test vector, non-symmetric form
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
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
@@ -2387,12 +2300,13 @@
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
endif
- enddo
- enddo
- enddo
- else
- stop 'wrong PML flag in PML memory variable calculation routine'
- endif
+ else
+ stop 'wrong PML flag in PML memory variable calculation routine'
+ endif
+ enddo
+ enddo
+ enddo
+
end subroutine pml_compute_memory_variables
More information about the CIG-COMMITS
mailing list