[cig-commits] r21035 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
joseph.charles at geodynamics.org
joseph.charles at geodynamics.org
Thu Nov 15 01:20:21 PST 2012
Author: joseph.charles
Date: 2012-11-15 01:20:21 -0800 (Thu, 15 Nov 2012)
New Revision: 21035
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_Dev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_Dev2.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90
Log:
fixes attenuation bug in CPU code by using 1st order Taylor expansion for epsilondev_loc variable
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90 2012-11-14 21:46:27 UTC (rev 21034)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90 2012-11-15 09:20:21 UTC (rev 21035)
@@ -61,13 +61,13 @@
else
! no optimizations used
- call compute_forces_elastic_noDev( iphase, NSPEC_AB,NGLOB_AB,displ,accel, &
+ call compute_forces_elastic_noDev( iphase, NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz, &
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,&
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION,&
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB, &
@@ -90,13 +90,13 @@
! adjoint simulations: backward/reconstructed wavefield
if( SIMULATION_TYPE == 3 ) &
call compute_forces_elastic_noDev( iphase, NSPEC_AB,NGLOB_AB,&
- b_displ,b_accel, &
+ b_displ,b_veloc,b_accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz, &
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,&
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION,&
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
b_alphaval,b_betaval,b_gammaval,&
NSPEC_ATTENUATION_AB, &
@@ -443,12 +443,12 @@
phase_ispec_inner_elastic,&
num_colors_outer_elastic,num_colors_inner_elastic)
#else
- call compute_forces_elastic_Dev_5p(iphase, NSPEC_AB,NGLOB_AB,displ,accel, &
+ call compute_forces_elastic_Dev_5p(iphase, NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
alphaval,betaval,gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -470,12 +470,12 @@
#endif
case (6)
- call compute_forces_elastic_Dev_6p(iphase, NSPEC_AB,NGLOB_AB,displ,accel, &
+ call compute_forces_elastic_Dev_6p(iphase, NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
alphaval,betaval,gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -496,12 +496,12 @@
phase_ispec_inner_elastic )
case (7)
- call compute_forces_elastic_Dev_7p(iphase, NSPEC_AB,NGLOB_AB,displ,accel, &
+ call compute_forces_elastic_Dev_7p(iphase, NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
alphaval,betaval,gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -522,12 +522,12 @@
phase_ispec_inner_elastic )
case (8)
- call compute_forces_elastic_Dev_8p(iphase, NSPEC_AB,NGLOB_AB,displ,accel, &
+ call compute_forces_elastic_Dev_8p(iphase, NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
alphaval,betaval,gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -548,12 +548,12 @@
phase_ispec_inner_elastic )
case (9)
- call compute_forces_elastic_Dev_9p(iphase, NSPEC_AB,NGLOB_AB,displ,accel, &
+ call compute_forces_elastic_Dev_9p(iphase, NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
alphaval,betaval,gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -574,12 +574,12 @@
phase_ispec_inner_elastic )
case (10)
- call compute_forces_elastic_Dev_10p(iphase, NSPEC_AB,NGLOB_AB,displ,accel, &
+ call compute_forces_elastic_Dev_10p(iphase, NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
alphaval,betaval,gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -629,12 +629,12 @@
case (5)
call compute_forces_elastic_Dev_5p(iphase, NSPEC_AB,NGLOB_AB, &
- b_displ,b_accel, &
+ b_displ,b_veloc,b_accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
b_alphaval,b_betaval,b_gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -656,12 +656,12 @@
case (6)
call compute_forces_elastic_Dev_6p(iphase, NSPEC_AB,NGLOB_AB, &
- b_displ,b_accel, &
+ b_displ,b_veloc,b_accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
b_alphaval,b_betaval,b_gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -683,12 +683,12 @@
case (7)
call compute_forces_elastic_Dev_7p(iphase, NSPEC_AB,NGLOB_AB, &
- b_displ,b_accel, &
+ b_displ,b_veloc,b_accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
b_alphaval,b_betaval,b_gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -710,12 +710,12 @@
case (8)
call compute_forces_elastic_Dev_8p(iphase, NSPEC_AB,NGLOB_AB, &
- b_displ,b_accel, &
+ b_displ,b_veloc,b_accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
b_alphaval,b_betaval,b_gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -737,12 +737,12 @@
case (9)
call compute_forces_elastic_Dev_9p(iphase, NSPEC_AB,NGLOB_AB, &
- b_displ,b_accel, &
+ b_displ,b_veloc,b_accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
b_alphaval,b_betaval,b_gammaval, &
NSPEC_ATTENUATION_AB, &
@@ -764,12 +764,12 @@
case (10)
call compute_forces_elastic_Dev_10p(iphase, NSPEC_AB,NGLOB_AB, &
- b_displ,b_accel, &
+ b_displ,b_veloc,b_accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
b_alphaval,b_betaval,b_gammaval, &
NSPEC_ATTENUATION_AB, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_Dev.f90 2012-11-14 21:46:27 UTC (rev 21034)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_Dev.f90 2012-11-15 09:20:21 UTC (rev 21035)
@@ -27,13 +27,13 @@
! Deville routine for NGLL == 5 (default)
subroutine compute_forces_elastic_Dev_5p( iphase ,NSPEC_AB,NGLOB_AB, &
- displ,accel, &
+ displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
@@ -62,9 +62,12 @@
integer :: NSPEC_AB,NGLOB_AB
-! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+! displacement, velocity and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
+! time step
+ real(kind=CUSTOM_REAL) :: deltat
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -124,6 +127,10 @@
real(kind=CUSTOM_REAL), dimension(5,5,5) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
! manually inline the calls to the Deville et al. (2002) routines
real(kind=CUSTOM_REAL), dimension(5,25) :: B1_m1_m2_5points,B2_m1_m2_5points,B3_m1_m2_5points
real(kind=CUSTOM_REAL), dimension(5,25) :: C1_m1_m2_5points,C2_m1_m2_5points,C3_m1_m2_5points
@@ -139,6 +146,19 @@
equivalence(newtempy1,E2_m1_m2_5points)
equivalence(newtempz1,E3_m1_m2_5points)
+ real(kind=CUSTOM_REAL), dimension(5,5,5) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+ real(kind=CUSTOM_REAL), dimension(5,5,5) :: dummyx_loc_att,dummyy_loc_att,dummyz_loc_att
+ real(kind=CUSTOM_REAL), dimension(5,25) :: B1_m1_m2_5points_att,B2_m1_m2_5points_att,B3_m1_m2_5points_att
+ real(kind=CUSTOM_REAL), dimension(5,25) :: C1_m1_m2_5points_att,C2_m1_m2_5points_att,C3_m1_m2_5points_att
+
+ equivalence(dummyx_loc_att,B1_m1_m2_5points_att)
+ equivalence(dummyy_loc_att,B2_m1_m2_5points_att)
+ equivalence(dummyz_loc_att,B3_m1_m2_5points_att)
+ equivalence(tempx1_att,C1_m1_m2_5points_att)
+ equivalence(tempy1_att,C2_m1_m2_5points_att)
+ equivalence(tempz1_att,C3_m1_m2_5points_att)
+
real(kind=CUSTOM_REAL), dimension(25,5) :: &
A1_mxm_m2_m1_5points,A2_mxm_m2_m1_5points,A3_mxm_m2_m1_5points
real(kind=CUSTOM_REAL), dimension(25,5) :: &
@@ -156,6 +176,18 @@
equivalence(newtempy3,E2_mxm_m2_m1_5points)
equivalence(newtempz3,E3_mxm_m2_m1_5points)
+ real(kind=CUSTOM_REAL), dimension(25,5) :: &
+ A1_mxm_m2_m1_5points_att,A2_mxm_m2_m1_5points_att,A3_mxm_m2_m1_5points_att
+ real(kind=CUSTOM_REAL), dimension(25,5) :: &
+ C1_mxm_m2_m1_5points_att,C2_mxm_m2_m1_5points_att,C3_mxm_m2_m1_5points_att
+
+ equivalence(dummyx_loc_att,A1_mxm_m2_m1_5points_att)
+ equivalence(dummyy_loc_att,A2_mxm_m2_m1_5points_att)
+ equivalence(dummyz_loc_att,A3_mxm_m2_m1_5points_att)
+ equivalence(tempx3_att,C1_mxm_m2_m1_5points_att)
+ equivalence(tempy3_att,C2_mxm_m2_m1_5points_att)
+ equivalence(tempz3_att,C3_mxm_m2_m1_5points_att)
+
! local attenuation parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
@@ -210,16 +242,31 @@
! stores displacment values in local array
do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
- dummyx_loc(i,j,k) = displ(1,iglob)
- dummyy_loc(i,j,k) = displ(2,iglob)
- dummyz_loc(i,j,k) = displ(3,iglob)
- enddo
- enddo
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc(i,j,k) = displ(1,iglob)
+ dummyy_loc(i,j,k) = displ(2,iglob)
+ dummyz_loc(i,j,k) = displ(3,iglob)
+ enddo
+ enddo
enddo
+ ! use first order Taylor expansion of displacement for local storage of stresses
+ ! at this current time step, to fix attenuation in a consistent way
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc_att(i,j,k) = deltat*veloc(1,iglob)
+ dummyy_loc_att(i,j,k) = deltat*veloc(2,iglob)
+ dummyz_loc_att(i,j,k) = deltat*veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+ endif
+
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
@@ -244,6 +291,34 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_5points_att(i,j) = C1_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B1_m1_m2_5points_att(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_5points_att(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_5points_att(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_5points_att(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_5points_att(5,j)
+
+ C2_m1_m2_5points_att(i,j) = C2_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B2_m1_m2_5points_att(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_5points_att(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_5points_att(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_5points_att(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_5points_att(5,j)
+
+ C3_m1_m2_5points_att(i,j) = C3_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B3_m1_m2_5points_att(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_5points_att(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_5points_att(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_5points_att(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_5points_att(5,j)
+ enddo
+ enddo
+ endif
+
! call mxm_m1_m1_5points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
! hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
do j=1,m1
@@ -269,6 +344,37 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2_att(i,j,k) = tempx2(i,j,k) + &
+ dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc_att(i,5,k)*hprime_xxT(5,j)
+
+ tempy2_att(i,j,k) = tempy2(i,j,k) + &
+ dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc_att(i,5,k)*hprime_xxT(5,j)
+
+ tempz2_att(i,j,k) = tempz2(i,j,k) + &
+ dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc_att(i,5,k)*hprime_xxT(5,j)
+ enddo
+ enddo
+ enddo
+ endif
+
! call mxm_m2_m1_5points(dummyx_loc,dummyy_loc,dummyz_loc,tempx3,tempy3,tempz3)
do j=1,m1
do i=1,m2
@@ -290,6 +396,34 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_5points_att(i,j) = C1_mxm_m2_m1_5points(i,j) + &
+ A1_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
+
+ C2_mxm_m2_m1_5points_att(i,j) = C2_mxm_m2_m1_5points(i,j) + &
+ A2_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
+
+ C3_mxm_m2_m1_5points_att(i,j) = C3_mxm_m2_m1_5points(i,j) + &
+ A3_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
+ enddo
+ enddo
+ endif
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -350,15 +484,44 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! computes deviatoric strain attenuation and/or for kernel calculations
- if (COMPUTE_AND_STORE_STRAIN) then
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
- epsilondev_xx_loc(i,j,k) = duxdxl - templ
- epsilondev_yy_loc(i,j,k) = duydyl - templ
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if ( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
+
+ ! compute deviatoric strain
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl_att - templ
+ epsilondev_yy_loc(i,j,k) = duydyl_att - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
kappal = kappastore(i,j,k,ispec)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_Dev2.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_Dev2.f90 2012-11-14 21:46:27 UTC (rev 21034)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_Dev2.f90 2012-11-15 09:20:21 UTC (rev 21035)
@@ -31,13 +31,13 @@
! for vectorizations when compiling
subroutine compute_forces_elastic_Dev_6p( iphase ,NSPEC_AB,NGLOB_AB, &
- displ,accel, &
+ displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
@@ -67,8 +67,11 @@
integer :: NSPEC_AB,NGLOB_AB
! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
+! time step
+ real(kind=CUSTOM_REAL) :: deltat
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -143,6 +146,19 @@
equivalence(newtempy1,E2_m1_m2_6points)
equivalence(newtempz1,E3_m1_m2_6points)
+ real(kind=CUSTOM_REAL), dimension(6,6,6) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+ real(kind=CUSTOM_REAL), dimension(6,6,6) :: dummyx_loc_att,dummyy_loc_att,dummyz_loc_att
+ real(kind=CUSTOM_REAL), dimension(6,36) :: B1_m1_m2_6points_att,B2_m1_m2_6points_att,B3_m1_m2_6points_att
+ real(kind=CUSTOM_REAL), dimension(6,36) :: C1_m1_m2_6points_att,C2_m1_m2_6points_att,C3_m1_m2_6points_att
+
+ equivalence(dummyx_loc_att,B1_m1_m2_6points_att)
+ equivalence(dummyy_loc_att,B2_m1_m2_6points_att)
+ equivalence(dummyz_loc_att,B3_m1_m2_6points_att)
+ equivalence(tempx1_att,C1_m1_m2_6points_att)
+ equivalence(tempy1_att,C2_m1_m2_6points_att)
+ equivalence(tempz1_att,C3_m1_m2_6points_att)
+
real(kind=CUSTOM_REAL), dimension(36,6) :: &
A1_mxm_m2_m1_6points,A2_mxm_m2_m1_6points,A3_mxm_m2_m1_6points
real(kind=CUSTOM_REAL), dimension(36,6) :: &
@@ -160,6 +176,18 @@
equivalence(newtempy3,E2_mxm_m2_m1_6points)
equivalence(newtempz3,E3_mxm_m2_m1_6points)
+ real(kind=CUSTOM_REAL), dimension(36,6) :: &
+ A1_mxm_m2_m1_6points_att,A2_mxm_m2_m1_6points_att,A3_mxm_m2_m1_6points_att
+ real(kind=CUSTOM_REAL), dimension(36,6) :: &
+ C1_mxm_m2_m1_6points_att,C2_mxm_m2_m1_6points_att,C3_mxm_m2_m1_6points_att
+
+ equivalence(dummyx_loc_att,A1_mxm_m2_m1_6points_att)
+ equivalence(dummyy_loc_att,A2_mxm_m2_m1_6points_att)
+ equivalence(dummyz_loc_att,A3_mxm_m2_m1_6points_att)
+ equivalence(tempx3_att,C1_mxm_m2_m1_6points_att)
+ equivalence(tempy3_att,C2_mxm_m2_m1_6points_att)
+ equivalence(tempz3_att,C3_mxm_m2_m1_6points_att)
+
! local attenuation parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
@@ -181,6 +209,10 @@
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
! local anisotropy parameters
real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
@@ -223,6 +255,21 @@
enddo
enddo
enddo
+
+ ! use first order Taylor expansion of displacement for local storage of stresses
+ ! at this current time step, to fix attenuation in a consistent way
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc_att(i,j,k) = deltat*veloc(1,iglob)
+ dummyy_loc_att(i,j,k) = deltat*veloc(2,iglob)
+ dummyz_loc_att(i,j,k) = deltat*veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+ endif
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
@@ -251,6 +298,37 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_6points_att(i,j) = C1_m1_m2_6points(i,j) + &
+ hprime_xx(i,1)*B1_m1_m2_6points_att(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_6points_att(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_6points_att(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_6points_att(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_6points_att(5,j) + &
+ hprime_xx(i,6)*B1_m1_m2_6points_att(6,j)
+
+ C2_m1_m2_6points_att(i,j) = C2_m1_m2_6points(i,j) + &
+ hprime_xx(i,1)*B2_m1_m2_6points_att(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_6points_att(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_6points_att(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_6points_att(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_6points_att(5,j) + &
+ hprime_xx(i,6)*B2_m1_m2_6points_att(6,j)
+
+ C3_m1_m2_6points_att(i,j) = C3_m1_m2_6points(i,j) + &
+ hprime_xx(i,1)*B3_m1_m2_6points_att(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_6points_att(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_6points_att(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_6points_att(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_6points_att(5,j) + &
+ hprime_xx(i,6)*B3_m1_m2_6points_att(6,j)
+ enddo
+ enddo
+ endif
+
! call mxm_m1_m1_6points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
! hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
do j=1,m1
@@ -279,6 +357,40 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2_att(i,j,k) = tempx2(i,j,k) + &
+ dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyx_loc_att(i,6,k)*hprime_xxT(6,j)
+
+ tempy2_att(i,j,k) = tempy2(i,j,k) + &
+ dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyy_loc_att(i,6,k)*hprime_xxT(6,j)
+
+ tempz2_att(i,j,k) = tempz2(i,j,k) + &
+ dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyz_loc_att(i,6,k)*hprime_xxT(6,j)
+ enddo
+ enddo
+ enddo
+ endif
+
! call mxm_m2_m1_6points(dummyx_loc,dummyy_loc,dummyz_loc,tempx3,tempy3,tempz3)
do j=1,m1
do i=1,m2
@@ -303,6 +415,37 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_6points_att(i,j) = C1_mxm_m2_m1_6points(i,j) + &
+ A1_mxm_m2_m1_6points_att(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_6points_att(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_6points_att(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_6points_att(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_6points_att(i,5)*hprime_xxT(5,j) + &
+ A1_mxm_m2_m1_6points_att(i,6)*hprime_xxT(6,j)
+
+ C2_mxm_m2_m1_6points_att(i,j) = C2_mxm_m2_m1_6points(i,j) + &
+ A2_mxm_m2_m1_6points_att(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_6points_att(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_6points_att(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_6points_att(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_6points_att(i,5)*hprime_xxT(5,j) + &
+ A2_mxm_m2_m1_6points_att(i,6)*hprime_xxT(6,j)
+
+ C3_mxm_m2_m1_6points_att(i,j) = C3_mxm_m2_m1_6points(i,j) + &
+ A3_mxm_m2_m1_6points_att(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_6points_att(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_6points_att(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_6points_att(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_6points_att(i,5)*hprime_xxT(5,j) + &
+ A3_mxm_m2_m1_6points_att(i,6)*hprime_xxT(6,j)
+ enddo
+ enddo
+ endif
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -363,15 +506,44 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! computes deviatoric strain attenuation and/or for kernel calculations
- if (COMPUTE_AND_STORE_STRAIN) then
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
- epsilondev_xx_loc(i,j,k) = duxdxl - templ
- epsilondev_yy_loc(i,j,k) = duydyl - templ
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if ( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
+
+ ! compute deviatoric strain
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl_att - templ
+ epsilondev_yy_loc(i,j,k) = duydyl_att - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
kappal = kappastore(i,j,k,ispec)
@@ -685,13 +857,13 @@
!
subroutine compute_forces_elastic_Dev_7p( iphase ,NSPEC_AB,NGLOB_AB, &
- displ,accel, &
+ displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
@@ -720,9 +892,12 @@
integer :: NSPEC_AB,NGLOB_AB
-! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+! displacement, velocity and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
+! time step
+ real(kind=CUSTOM_REAL) :: deltat
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -797,6 +972,19 @@
equivalence(newtempy1,E2_m1_m2_7points)
equivalence(newtempz1,E3_m1_m2_7points)
+ real(kind=CUSTOM_REAL), dimension(7,7,7) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+ real(kind=CUSTOM_REAL), dimension(7,7,7) :: dummyx_loc_att,dummyy_loc_att,dummyz_loc_att
+ real(kind=CUSTOM_REAL), dimension(7,49) :: B1_m1_m2_7points_att,B2_m1_m2_7points_att,B3_m1_m2_7points_att
+ real(kind=CUSTOM_REAL), dimension(7,49) :: C1_m1_m2_7points_att,C2_m1_m2_7points_att,C3_m1_m2_7points_att
+
+ equivalence(dummyx_loc_att,B1_m1_m2_7points_att)
+ equivalence(dummyy_loc_att,B2_m1_m2_7points_att)
+ equivalence(dummyz_loc_att,B3_m1_m2_7points_att)
+ equivalence(tempx1_att,C1_m1_m2_7points_att)
+ equivalence(tempy1_att,C2_m1_m2_7points_att)
+ equivalence(tempz1_att,C3_m1_m2_7points_att)
+
real(kind=CUSTOM_REAL), dimension(49,7) :: &
A1_mxm_m2_m1_7points,A2_mxm_m2_m1_7points,A3_mxm_m2_m1_7points
real(kind=CUSTOM_REAL), dimension(49,7) :: &
@@ -814,6 +1002,18 @@
equivalence(newtempy3,E2_mxm_m2_m1_7points)
equivalence(newtempz3,E3_mxm_m2_m1_7points)
+ real(kind=CUSTOM_REAL), dimension(49,7) :: &
+ A1_mxm_m2_m1_7points_att,A2_mxm_m2_m1_7points_att,A3_mxm_m2_m1_7points_att
+ real(kind=CUSTOM_REAL), dimension(49,7) :: &
+ C1_mxm_m2_m1_7points_att,C2_mxm_m2_m1_7points_att,C3_mxm_m2_m1_7points_att
+
+ equivalence(dummyx_loc_att,A1_mxm_m2_m1_7points_att)
+ equivalence(dummyy_loc_att,A2_mxm_m2_m1_7points_att)
+ equivalence(dummyz_loc_att,A3_mxm_m2_m1_7points_att)
+ equivalence(tempx3_att,C1_mxm_m2_m1_7points_att)
+ equivalence(tempy3_att,C2_mxm_m2_m1_7points_att)
+ equivalence(tempz3_att,C3_mxm_m2_m1_7points_att)
+
! local attenuation parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
@@ -835,6 +1035,10 @@
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
! local anisotropy parameters
real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
@@ -878,6 +1082,21 @@
enddo
enddo
+ ! use first order Taylor expansion of displacement for local storage of stresses
+ ! at this current time step, to fix attenuation in a consistent way
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc_att(i,j,k) = deltat*veloc(1,iglob)
+ dummyy_loc_att(i,j,k) = deltat*veloc(2,iglob)
+ dummyz_loc_att(i,j,k) = deltat*veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+ endif
+
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
@@ -908,6 +1127,40 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_7points_att(i,j) = C1_m1_m2_7points(i,j) + &
+ hprime_xx(i,1)*B1_m1_m2_7points_att(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_7points_att(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_7points_att(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_7points_att(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_7points_att(5,j) + &
+ hprime_xx(i,6)*B1_m1_m2_7points_att(6,j) + &
+ hprime_xx(i,7)*B1_m1_m2_7points_att(7,j)
+
+ C2_m1_m2_7points_att(i,j) = C2_m1_m2_7points(i,j) + &
+ hprime_xx(i,1)*B2_m1_m2_7points_att(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_7points_att(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_7points_att(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_7points_att(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_7points_att(5,j) + &
+ hprime_xx(i,6)*B2_m1_m2_7points_att(6,j) + &
+ hprime_xx(i,7)*B2_m1_m2_7points_att(7,j)
+
+ C3_m1_m2_7points_att(i,j) = C3_m1_m2_7points(i,j) + &
+ hprime_xx(i,1)*B3_m1_m2_7points_att(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_7points_att(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_7points_att(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_7points_att(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_7points_att(5,j) + &
+ hprime_xx(i,6)*B3_m1_m2_7points_att(6,j) + &
+ hprime_xx(i,7)*B3_m1_m2_7points_att(7,j)
+ enddo
+ enddo
+ endif
+
! call mxm_m1_m1_7points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
! hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
do j=1,m1
@@ -939,6 +1192,43 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2_att(i,j,k) = tempx2(i,j,k) + &
+ dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyx_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyx_loc_att(i,7,k)*hprime_xxT(7,j)
+
+ tempy2_att(i,j,k) = tempy2(i,j,k) + &
+ dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyy_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyy_loc_att(i,7,k)*hprime_xxT(7,j)
+
+ tempz2_att(i,j,k) = tempz2(i,j,k) + &
+ dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyz_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyz_loc_att(i,7,k)*hprime_xxT(7,j)
+ enddo
+ enddo
+ enddo
+ endif
+
! call mxm_m2_m1_7points(dummyx_loc,dummyy_loc,dummyz_loc,tempx3,tempy3,tempz3)
do j=1,m1
do i=1,m2
@@ -966,6 +1256,40 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_7points_att(i,j) = C1_mxm_m2_m1_7points(i,j) + &
+ A1_mxm_m2_m1_7points_att(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_7points_att(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_7points_att(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_7points_att(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_7points_att(i,5)*hprime_xxT(5,j) + &
+ A1_mxm_m2_m1_7points_att(i,6)*hprime_xxT(6,j) + &
+ A1_mxm_m2_m1_7points_att(i,7)*hprime_xxT(7,j)
+
+ C2_mxm_m2_m1_7points_att(i,j) = C2_mxm_m2_m1_7points(i,j) + &
+ A2_mxm_m2_m1_7points_att(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_7points_att(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_7points_att(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_7points_att(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_7points_att(i,5)*hprime_xxT(5,j) + &
+ A2_mxm_m2_m1_7points_att(i,6)*hprime_xxT(6,j) + &
+ A2_mxm_m2_m1_7points_att(i,7)*hprime_xxT(7,j)
+
+ C3_mxm_m2_m1_7points_att(i,j) = C3_mxm_m2_m1_7points(i,j) + &
+ A3_mxm_m2_m1_7points_att(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_7points_att(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_7points_att(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_7points_att(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_7points_att(i,5)*hprime_xxT(5,j) + &
+ A3_mxm_m2_m1_7points_att(i,6)*hprime_xxT(6,j) + &
+ A3_mxm_m2_m1_7points_att(i,7)*hprime_xxT(7,j)
+ enddo
+ enddo
+ endif
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -1026,15 +1350,44 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! computes deviatoric strain attenuation and/or for kernel calculations
- if (COMPUTE_AND_STORE_STRAIN) then
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
- epsilondev_xx_loc(i,j,k) = duxdxl - templ
- epsilondev_yy_loc(i,j,k) = duydyl - templ
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if ( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
+
+ ! compute deviatoric strain
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl_att - templ
+ epsilondev_yy_loc(i,j,k) = duydyl_att - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
kappal = kappastore(i,j,k,ispec)
@@ -1357,13 +1710,13 @@
!
subroutine compute_forces_elastic_Dev_8p( iphase ,NSPEC_AB,NGLOB_AB, &
- displ,accel, &
+ displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
@@ -1392,9 +1745,12 @@
integer :: NSPEC_AB,NGLOB_AB
-! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+! displacement, velocity and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
+! time step
+ real(kind=CUSTOM_REAL) :: deltat
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -1469,6 +1825,19 @@
equivalence(newtempy1,E2_m1_m2_8points)
equivalence(newtempz1,E3_m1_m2_8points)
+ real(kind=CUSTOM_REAL), dimension(8,8,8) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+ real(kind=CUSTOM_REAL), dimension(8,8,8) :: dummyx_loc_att,dummyy_loc_att,dummyz_loc_att
+ real(kind=CUSTOM_REAL), dimension(8,64) :: B1_m1_m2_8points_att,B2_m1_m2_8points_att,B3_m1_m2_8points_att
+ real(kind=CUSTOM_REAL), dimension(8,64) :: C1_m1_m2_8points_att,C2_m1_m2_8points_att,C3_m1_m2_8points_att
+
+ equivalence(dummyx_loc_att,B1_m1_m2_8points_att)
+ equivalence(dummyy_loc_att,B2_m1_m2_8points_att)
+ equivalence(dummyz_loc_att,B3_m1_m2_8points_att)
+ equivalence(tempx1_att,C1_m1_m2_8points_att)
+ equivalence(tempy1_att,C2_m1_m2_8points_att)
+ equivalence(tempz1_att,C3_m1_m2_8points_att)
+
real(kind=CUSTOM_REAL), dimension(64,8) :: &
A1_mxm_m2_m1_8points,A2_mxm_m2_m1_8points,A3_mxm_m2_m1_8points
real(kind=CUSTOM_REAL), dimension(64,8) :: &
@@ -1486,6 +1855,18 @@
equivalence(newtempy3,E2_mxm_m2_m1_8points)
equivalence(newtempz3,E3_mxm_m2_m1_8points)
+ real(kind=CUSTOM_REAL), dimension(64,8) :: &
+ A1_mxm_m2_m1_8points_att,A2_mxm_m2_m1_8points_att,A3_mxm_m2_m1_8points_att
+ real(kind=CUSTOM_REAL), dimension(64,8) :: &
+ C1_mxm_m2_m1_8points_att,C2_mxm_m2_m1_8points_att,C3_mxm_m2_m1_8points_att
+
+ equivalence(dummyx_loc_att,A1_mxm_m2_m1_8points_att)
+ equivalence(dummyy_loc_att,A2_mxm_m2_m1_8points_att)
+ equivalence(dummyz_loc_att,A3_mxm_m2_m1_8points_att)
+ equivalence(tempx3_att,C1_mxm_m2_m1_8points_att)
+ equivalence(tempy3_att,C2_mxm_m2_m1_8points_att)
+ equivalence(tempz3_att,C3_mxm_m2_m1_8points_att)
+
! local attenuation parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
@@ -1507,6 +1888,10 @@
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
! local anisotropy parameters
real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
@@ -1550,6 +1935,21 @@
enddo
enddo
+ ! use first order Taylor expansion of displacement for local storage of stresses
+ ! at this current time step, to fix attenuation in a consistent way
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc_att(i,j,k) = deltat*veloc(1,iglob)
+ dummyy_loc_att(i,j,k) = deltat*veloc(2,iglob)
+ dummyz_loc_att(i,j,k) = deltat*veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+ endif
+
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
@@ -1583,6 +1983,43 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_8points_att(i,j) = C1_m1_m2_8points(i,j) + &
+ hprime_xx(i,1)*B1_m1_m2_8points_att(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_8points_att(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_8points_att(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_8points_att(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_8points_att(5,j) + &
+ hprime_xx(i,6)*B1_m1_m2_8points_att(6,j) + &
+ hprime_xx(i,7)*B1_m1_m2_8points_att(7,j) + &
+ hprime_xx(i,8)*B1_m1_m2_8points_att(8,j)
+
+ C2_m1_m2_8points_att(i,j) = C2_m1_m2_8points(i,j) + &
+ hprime_xx(i,1)*B2_m1_m2_8points_att(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_8points_att(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_8points_att(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_8points_att(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_8points_att(5,j) + &
+ hprime_xx(i,6)*B2_m1_m2_8points_att(6,j) + &
+ hprime_xx(i,7)*B2_m1_m2_8points_att(7,j) + &
+ hprime_xx(i,8)*B2_m1_m2_8points_att(8,j)
+
+ C3_m1_m2_8points_att(i,j) = C3_m1_m2_8points(i,j) + &
+ hprime_xx(i,1)*B3_m1_m2_8points_att(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_8points_att(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_8points_att(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_8points_att(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_8points_att(5,j) + &
+ hprime_xx(i,6)*B3_m1_m2_8points_att(6,j) + &
+ hprime_xx(i,7)*B3_m1_m2_8points_att(7,j) + &
+ hprime_xx(i,8)*B3_m1_m2_8points_att(8,j)
+ enddo
+ enddo
+ endif
+
! call mxm_m1_m1_8points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
! hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
do j=1,m1
@@ -1617,6 +2054,46 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2_att(i,j,k) = tempx2(i,j,k) + &
+ dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyx_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyx_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyx_loc_att(i,8,k)*hprime_xxT(8,j)
+
+ tempy2_att(i,j,k) = tempy2(i,j,k) + &
+ dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyy_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyy_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyy_loc_att(i,8,k)*hprime_xxT(8,j)
+
+ tempz2_att(i,j,k) = tempz2(i,j,k) + &
+ dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyz_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyz_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyz_loc_att(i,8,k)*hprime_xxT(8,j)
+ enddo
+ enddo
+ enddo
+ endif
+
! call mxm_m2_m1_8points(dummyx_loc,dummyy_loc,dummyz_loc,tempx3,tempy3,tempz3)
do j=1,m1
do i=1,m2
@@ -1647,6 +2124,43 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_8points_att(i,j) = C1_mxm_m2_m1_8points(i,j) + &
+ A1_mxm_m2_m1_8points_att(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_8points_att(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_8points_att(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_8points_att(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_8points_att(i,5)*hprime_xxT(5,j) + &
+ A1_mxm_m2_m1_8points_att(i,6)*hprime_xxT(6,j) + &
+ A1_mxm_m2_m1_8points_att(i,7)*hprime_xxT(7,j) + &
+ A1_mxm_m2_m1_8points_att(i,8)*hprime_xxT(8,j)
+
+ C2_mxm_m2_m1_8points_att(i,j) = C2_mxm_m2_m1_8points(i,j) + &
+ A2_mxm_m2_m1_8points_att(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_8points_att(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_8points_att(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_8points_att(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_8points_att(i,5)*hprime_xxT(5,j) + &
+ A2_mxm_m2_m1_8points_att(i,6)*hprime_xxT(6,j) + &
+ A2_mxm_m2_m1_8points_att(i,7)*hprime_xxT(7,j) + &
+ A2_mxm_m2_m1_8points_att(i,8)*hprime_xxT(8,j)
+
+ C3_mxm_m2_m1_8points_att(i,j) = C3_mxm_m2_m1_8points(i,j) + &
+ A3_mxm_m2_m1_8points_att(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_8points_att(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_8points_att(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_8points_att(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_8points_att(i,5)*hprime_xxT(5,j) + &
+ A3_mxm_m2_m1_8points_att(i,6)*hprime_xxT(6,j) + &
+ A3_mxm_m2_m1_8points_att(i,7)*hprime_xxT(7,j) + &
+ A3_mxm_m2_m1_8points_att(i,8)*hprime_xxT(8,j)
+ enddo
+ enddo
+ endif
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -1707,15 +2221,44 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! computes deviatoric strain attenuation and/or for kernel calculations
- if (COMPUTE_AND_STORE_STRAIN) then
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
- epsilondev_xx_loc(i,j,k) = duxdxl - templ
- epsilondev_yy_loc(i,j,k) = duydyl - templ
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if ( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
+
+ ! compute deviatoric strain
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl_att - templ
+ epsilondev_yy_loc(i,j,k) = duydyl_att - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
kappal = kappastore(i,j,k,ispec)
@@ -2047,13 +2590,13 @@
!
subroutine compute_forces_elastic_Dev_9p( iphase ,NSPEC_AB,NGLOB_AB, &
- displ,accel, &
+ displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
@@ -2082,9 +2625,12 @@
integer :: NSPEC_AB,NGLOB_AB
-! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+! displacement, velocity and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
+! time step
+ real(kind=CUSTOM_REAL) :: deltat
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -2159,6 +2705,19 @@
equivalence(newtempy1,E2_m1_m2_9points)
equivalence(newtempz1,E3_m1_m2_9points)
+ real(kind=CUSTOM_REAL), dimension(9,9,9) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+ real(kind=CUSTOM_REAL), dimension(9,9,9) :: dummyx_loc_att,dummyy_loc_att,dummyz_loc_att
+ real(kind=CUSTOM_REAL), dimension(9,81) :: B1_m1_m2_9points_att,B2_m1_m2_9points_att,B3_m1_m2_9points_att
+ real(kind=CUSTOM_REAL), dimension(9,81) :: C1_m1_m2_9points_att,C2_m1_m2_9points_att,C3_m1_m2_9points_att
+
+ equivalence(dummyx_loc_att,B1_m1_m2_9points_att)
+ equivalence(dummyy_loc_att,B2_m1_m2_9points_att)
+ equivalence(dummyz_loc_att,B3_m1_m2_9points_att)
+ equivalence(tempx1_att,C1_m1_m2_9points_att)
+ equivalence(tempy1_att,C2_m1_m2_9points_att)
+ equivalence(tempz1_att,C3_m1_m2_9points_att)
+
real(kind=CUSTOM_REAL), dimension(81,9) :: &
A1_mxm_m2_m1_9points,A2_mxm_m2_m1_9points,A3_mxm_m2_m1_9points
real(kind=CUSTOM_REAL), dimension(81,9) :: &
@@ -2176,6 +2735,18 @@
equivalence(newtempy3,E2_mxm_m2_m1_9points)
equivalence(newtempz3,E3_mxm_m2_m1_9points)
+ real(kind=CUSTOM_REAL), dimension(81,9) :: &
+ A1_mxm_m2_m1_9points_att,A2_mxm_m2_m1_9points_att,A3_mxm_m2_m1_9points_att
+ real(kind=CUSTOM_REAL), dimension(81,9) :: &
+ C1_mxm_m2_m1_9points_att,C2_mxm_m2_m1_9points_att,C3_mxm_m2_m1_9points_att
+
+ equivalence(dummyx_loc_att,A1_mxm_m2_m1_9points_att)
+ equivalence(dummyy_loc_att,A2_mxm_m2_m1_9points_att)
+ equivalence(dummyz_loc_att,A3_mxm_m2_m1_9points_att)
+ equivalence(tempx3_att,C1_mxm_m2_m1_9points_att)
+ equivalence(tempy3_att,C2_mxm_m2_m1_9points_att)
+ equivalence(tempz3_att,C3_mxm_m2_m1_9points_att)
+
! local attenuation parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
@@ -2197,6 +2768,10 @@
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
! local anisotropy parameters
real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
@@ -2240,6 +2815,21 @@
enddo
enddo
+ ! use first order Taylor expansion of displacement for local storage of stresses
+ ! at this current time step, to fix attenuation in a consistent way
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc_att(i,j,k) = deltat*veloc(1,iglob)
+ dummyy_loc_att(i,j,k) = deltat*veloc(2,iglob)
+ dummyz_loc_att(i,j,k) = deltat*veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+ endif
+
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
@@ -2276,6 +2866,46 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_9points_att(i,j) = C1_m1_m2_9points(i,j) + &
+ hprime_xx(i,1)*B1_m1_m2_9points_att(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_9points_att(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_9points_att(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_9points_att(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_9points_att(5,j) + &
+ hprime_xx(i,6)*B1_m1_m2_9points_att(6,j) + &
+ hprime_xx(i,7)*B1_m1_m2_9points_att(7,j) + &
+ hprime_xx(i,8)*B1_m1_m2_9points_att(8,j) + &
+ hprime_xx(i,9)*B1_m1_m2_9points_att(9,j)
+
+ C2_m1_m2_9points_att(i,j) = C2_m1_m2_9points(i,j) + &
+ hprime_xx(i,1)*B2_m1_m2_9points_att(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_9points_att(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_9points_att(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_9points_att(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_9points_att(5,j) + &
+ hprime_xx(i,6)*B2_m1_m2_9points_att(6,j) + &
+ hprime_xx(i,7)*B2_m1_m2_9points_att(7,j) + &
+ hprime_xx(i,8)*B2_m1_m2_9points_att(8,j) + &
+ hprime_xx(i,9)*B2_m1_m2_9points_att(9,j)
+
+ C3_m1_m2_9points_att(i,j) = C3_m1_m2_9points(i,j) + &
+ hprime_xx(i,1)*B3_m1_m2_9points_att(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_9points_att(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_9points_att(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_9points_att(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_9points_att(5,j) + &
+ hprime_xx(i,6)*B3_m1_m2_9points_att(6,j) + &
+ hprime_xx(i,7)*B3_m1_m2_9points_att(7,j) + &
+ hprime_xx(i,8)*B3_m1_m2_9points_att(8,j) + &
+ hprime_xx(i,9)*B3_m1_m2_9points_att(9,j)
+ enddo
+ enddo
+ endif
+
! call mxm_m1_m1_9points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
! hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
do j=1,m1
@@ -2313,6 +2943,49 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2_att(i,j,k) = tempx2(i,j,k) + &
+ dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyx_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyx_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyx_loc_att(i,8,k)*hprime_xxT(8,j) + &
+ dummyx_loc_att(i,9,k)*hprime_xxT(9,j)
+
+ tempy2_att(i,j,k) = tempy2(i,j,k) + &
+ dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyy_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyy_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyy_loc_att(i,8,k)*hprime_xxT(8,j) + &
+ dummyy_loc_att(i,9,k)*hprime_xxT(9,j)
+
+ tempz2_att(i,j,k) = tempz2(i,j,k) + &
+ dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyz_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyz_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyz_loc_att(i,8,k)*hprime_xxT(8,j) + &
+ dummyz_loc_att(i,9,k)*hprime_xxT(9,j)
+ enddo
+ enddo
+ enddo
+ endif
+
! call mxm_m2_m1_9points(dummyx_loc,dummyy_loc,dummyz_loc,tempx3,tempy3,tempz3)
do j=1,m1
do i=1,m2
@@ -2346,6 +3019,46 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_9points_att(i,j) = C1_mxm_m2_m1_9points(i,j) + &
+ A1_mxm_m2_m1_9points_att(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_9points_att(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_9points_att(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_9points_att(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_9points_att(i,5)*hprime_xxT(5,j) + &
+ A1_mxm_m2_m1_9points_att(i,6)*hprime_xxT(6,j) + &
+ A1_mxm_m2_m1_9points_att(i,7)*hprime_xxT(7,j) + &
+ A1_mxm_m2_m1_9points_att(i,8)*hprime_xxT(8,j) + &
+ A1_mxm_m2_m1_9points_att(i,9)*hprime_xxT(9,j)
+
+ C2_mxm_m2_m1_9points_att(i,j) = C2_mxm_m2_m1_9points(i,j) + &
+ A2_mxm_m2_m1_9points_att(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_9points_att(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_9points_att(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_9points_att(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_9points_att(i,5)*hprime_xxT(5,j) + &
+ A2_mxm_m2_m1_9points_att(i,6)*hprime_xxT(6,j) + &
+ A2_mxm_m2_m1_9points_att(i,7)*hprime_xxT(7,j) + &
+ A2_mxm_m2_m1_9points_att(i,8)*hprime_xxT(8,j) + &
+ A2_mxm_m2_m1_9points_att(i,9)*hprime_xxT(9,j)
+
+ C3_mxm_m2_m1_9points_att(i,j) = C3_mxm_m2_m1_9points(i,j) + &
+ A3_mxm_m2_m1_9points_att(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_9points_att(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_9points_att(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_9points_att(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_9points_att(i,5)*hprime_xxT(5,j) + &
+ A3_mxm_m2_m1_9points_att(i,6)*hprime_xxT(6,j) + &
+ A3_mxm_m2_m1_9points_att(i,7)*hprime_xxT(7,j) + &
+ A3_mxm_m2_m1_9points_att(i,8)*hprime_xxT(8,j) + &
+ A3_mxm_m2_m1_9points_att(i,9)*hprime_xxT(9,j)
+ enddo
+ enddo
+ endif
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -2406,15 +3119,44 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! computes deviatoric strain attenuation and/or for kernel calculations
- if (COMPUTE_AND_STORE_STRAIN) then
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
- epsilondev_xx_loc(i,j,k) = duxdxl - templ
- epsilondev_yy_loc(i,j,k) = duydyl - templ
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if ( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
+
+ ! compute deviatoric strain
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl_att - templ
+ epsilondev_yy_loc(i,j,k) = duydyl_att - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
kappal = kappastore(i,j,k,ispec)
@@ -2755,13 +3497,13 @@
!
subroutine compute_forces_elastic_Dev_10p( iphase ,NSPEC_AB,NGLOB_AB, &
- displ,accel, &
+ displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
@@ -2790,9 +3532,12 @@
integer :: NSPEC_AB,NGLOB_AB
-! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+! displacement, velocity and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
+! time step
+ real(kind=CUSTOM_REAL) :: deltat
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -2867,6 +3612,19 @@
equivalence(newtempy1,E2_m1_m2_10points)
equivalence(newtempz1,E3_m1_m2_10points)
+ real(kind=CUSTOM_REAL), dimension(10,10,10) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+ real(kind=CUSTOM_REAL), dimension(10,10,10) :: dummyx_loc_att,dummyy_loc_att,dummyz_loc_att
+ real(kind=CUSTOM_REAL), dimension(10,100) :: B1_m1_m2_10points_att,B2_m1_m2_10points_att,B3_m1_m2_10points_att
+ real(kind=CUSTOM_REAL), dimension(10,100) :: C1_m1_m2_10points_att,C2_m1_m2_10points_att,C3_m1_m2_10points_att
+
+ equivalence(dummyx_loc_att,B1_m1_m2_10points_att)
+ equivalence(dummyy_loc_att,B2_m1_m2_10points_att)
+ equivalence(dummyz_loc_att,B3_m1_m2_10points_att)
+ equivalence(tempx1_att,C1_m1_m2_10points_att)
+ equivalence(tempy1_att,C2_m1_m2_10points_att)
+ equivalence(tempz1_att,C3_m1_m2_10points_att)
+
real(kind=CUSTOM_REAL), dimension(100,10) :: &
A1_mxm_m2_m1_10points,A2_mxm_m2_m1_10points,A3_mxm_m2_m1_10points
real(kind=CUSTOM_REAL), dimension(100,10) :: &
@@ -2884,6 +3642,18 @@
equivalence(newtempy3,E2_mxm_m2_m1_10points)
equivalence(newtempz3,E3_mxm_m2_m1_10points)
+ real(kind=CUSTOM_REAL), dimension(100,10) :: &
+ A1_mxm_m2_m1_10points_att,A2_mxm_m2_m1_10points_att,A3_mxm_m2_m1_10points_att
+ real(kind=CUSTOM_REAL), dimension(100,10) :: &
+ C1_mxm_m2_m1_10points_att,C2_mxm_m2_m1_10points_att,C3_mxm_m2_m1_10points_att
+
+ equivalence(dummyx_loc_att,A1_mxm_m2_m1_10points_att)
+ equivalence(dummyy_loc_att,A2_mxm_m2_m1_10points_att)
+ equivalence(dummyz_loc_att,A3_mxm_m2_m1_10points_att)
+ equivalence(tempx3_att,C1_mxm_m2_m1_10points_att)
+ equivalence(tempy3_att,C2_mxm_m2_m1_10points_att)
+ equivalence(tempz3_att,C3_mxm_m2_m1_10points_att)
+
! local attenuation parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
@@ -2905,6 +3675,10 @@
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
! local anisotropy parameters
real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
@@ -2948,6 +3722,21 @@
enddo
enddo
+ ! use first order Taylor expansion of displacement for local storage of stresses
+ ! at this current time step, to fix attenuation in a consistent way
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc_att(i,j,k) = deltat*veloc(1,iglob)
+ dummyy_loc_att(i,j,k) = deltat*veloc(2,iglob)
+ dummyz_loc_att(i,j,k) = deltat*veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+ endif
+
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
@@ -2987,6 +3776,49 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_10points_att(i,j) = C1_m1_m2_10points(i,j) + &
+ hprime_xx(i,1)*B1_m1_m2_10points_att(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_10points_att(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_10points_att(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_10points_att(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_10points_att(5,j) + &
+ hprime_xx(i,6)*B1_m1_m2_10points_att(6,j) + &
+ hprime_xx(i,7)*B1_m1_m2_10points_att(7,j) + &
+ hprime_xx(i,8)*B1_m1_m2_10points_att(8,j) + &
+ hprime_xx(i,9)*B1_m1_m2_10points_att(9,j) + &
+ hprime_xx(i,10)*B1_m1_m2_10points_att(10,j)
+
+ C2_m1_m2_10points_att(i,j) = C2_m1_m2_10points(i,j) + &
+ hprime_xx(i,1)*B2_m1_m2_10points_att(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_10points_att(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_10points_att(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_10points_att(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_10points_att(5,j) + &
+ hprime_xx(i,6)*B2_m1_m2_10points_att(6,j) + &
+ hprime_xx(i,7)*B2_m1_m2_10points_att(7,j) + &
+ hprime_xx(i,8)*B2_m1_m2_10points_att(8,j) + &
+ hprime_xx(i,9)*B2_m1_m2_10points_att(9,j) + &
+ hprime_xx(i,10)*B2_m1_m2_10points_att(10,j)
+
+ C3_m1_m2_10points_att(i,j) = C3_m1_m2_10points(i,j) + &
+ hprime_xx(i,1)*B3_m1_m2_10points_att(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_10points_att(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_10points_att(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_10points_att(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_10points_att(5,j) + &
+ hprime_xx(i,6)*B3_m1_m2_10points_att(6,j) + &
+ hprime_xx(i,7)*B3_m1_m2_10points_att(7,j) + &
+ hprime_xx(i,8)*B3_m1_m2_10points_att(8,j) + &
+ hprime_xx(i,9)*B3_m1_m2_10points_att(9,j) + &
+ hprime_xx(i,10)*B3_m1_m2_10points_att(10,j)
+ enddo
+ enddo
+ endif
+
! call mxm_m1_m1_10points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
! hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
do j=1,m1
@@ -3027,6 +3859,52 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2_att(i,j,k) = tempx2(i,j,k) + &
+ dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyx_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyx_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyx_loc_att(i,8,k)*hprime_xxT(8,j) + &
+ dummyx_loc_att(i,9,k)*hprime_xxT(9,j) + &
+ dummyx_loc_att(i,10,k)*hprime_xxT(10,j)
+
+ tempy2_att(i,j,k) = tempy2(i,j,k) + &
+ dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyy_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyy_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyy_loc_att(i,8,k)*hprime_xxT(8,j) + &
+ dummyy_loc_att(i,9,k)*hprime_xxT(9,j) + &
+ dummyy_loc_att(i,10,k)*hprime_xxT(10,j)
+
+ tempz2_att(i,j,k) = tempz2(i,j,k) + &
+ dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc_att(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc_att(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc_att(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc_att(i,5,k)*hprime_xxT(5,j) + &
+ dummyz_loc_att(i,6,k)*hprime_xxT(6,j) + &
+ dummyz_loc_att(i,7,k)*hprime_xxT(7,j) + &
+ dummyz_loc_att(i,8,k)*hprime_xxT(8,j) + &
+ dummyz_loc_att(i,9,k)*hprime_xxT(9,j) + &
+ dummyz_loc_att(i,10,k)*hprime_xxT(10,j)
+ enddo
+ enddo
+ enddo
+ endif
+
! call mxm_m2_m1_10points(dummyx_loc,dummyy_loc,dummyz_loc,tempx3,tempy3,tempz3)
do j=1,m1
do i=1,m2
@@ -3063,6 +3941,49 @@
enddo
enddo
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_10points_att(i,j) = C1_mxm_m2_m1_10points(i,j) + &
+ A1_mxm_m2_m1_10points_att(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_10points_att(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_10points_att(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_10points_att(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_10points_att(i,5)*hprime_xxT(5,j) + &
+ A1_mxm_m2_m1_10points_att(i,6)*hprime_xxT(6,j) + &
+ A1_mxm_m2_m1_10points_att(i,7)*hprime_xxT(7,j) + &
+ A1_mxm_m2_m1_10points_att(i,8)*hprime_xxT(8,j) + &
+ A1_mxm_m2_m1_10points_att(i,9)*hprime_xxT(9,j) + &
+ A1_mxm_m2_m1_10points_att(i,10)*hprime_xxT(10,j)
+
+ C2_mxm_m2_m1_10points_att(i,j) = C2_mxm_m2_m1_10points(i,j) + &
+ A2_mxm_m2_m1_10points_att(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_10points_att(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_10points_att(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_10points_att(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_10points_att(i,5)*hprime_xxT(5,j) + &
+ A2_mxm_m2_m1_10points_att(i,6)*hprime_xxT(6,j) + &
+ A2_mxm_m2_m1_10points_att(i,7)*hprime_xxT(7,j) + &
+ A2_mxm_m2_m1_10points_att(i,8)*hprime_xxT(8,j) + &
+ A2_mxm_m2_m1_10points_att(i,9)*hprime_xxT(9,j) + &
+ A2_mxm_m2_m1_10points_att(i,10)*hprime_xxT(10,j)
+
+ C3_mxm_m2_m1_10points_att(i,j) = C3_mxm_m2_m1_10points(i,j) + &
+ A3_mxm_m2_m1_10points_att(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_10points_att(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_10points_att(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_10points_att(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_10points_att(i,5)*hprime_xxT(5,j) + &
+ A3_mxm_m2_m1_10points_att(i,6)*hprime_xxT(6,j) + &
+ A3_mxm_m2_m1_10points_att(i,7)*hprime_xxT(7,j) + &
+ A3_mxm_m2_m1_10points_att(i,8)*hprime_xxT(8,j) + &
+ A3_mxm_m2_m1_10points_att(i,9)*hprime_xxT(9,j) + &
+ A3_mxm_m2_m1_10points_att(i,10)*hprime_xxT(10,j)
+ enddo
+ enddo
+ endif
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -3123,15 +4044,45 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! computes deviatoric strain attenuation and/or for kernel calculations
- if (COMPUTE_AND_STORE_STRAIN) then
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
- epsilondev_xx_loc(i,j,k) = duxdxl - templ
- epsilondev_yy_loc(i,j,k) = duydyl - templ
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+
+ if ( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
+
+ ! compute deviatoric strain
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl_att - templ
+ epsilondev_yy_loc(i,j,k) = duydyl_att - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
kappal = kappastore(i,j,k,ispec)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90 2012-11-14 21:46:27 UTC (rev 21034)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90 2012-11-15 09:20:21 UTC (rev 21035)
@@ -25,30 +25,30 @@
!=====================================================================
subroutine compute_forces_elastic_noDev( iphase, &
- NSPEC_AB,NGLOB_AB,displ,accel,&
+ NSPEC_AB,NGLOB_AB,displ,veloc,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- hprime_xx,hprime_yy,hprime_zz,&
+ hprime_xx,hprime_yy,hprime_zz, &
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,&
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- kappastore,mustore,jacobian,ibool,&
- ATTENUATION,&
+ kappastore,mustore,jacobian,ibool, &
+ ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
- alphaval,betaval,gammaval,&
+ alphaval,betaval,gammaval, &
NSPEC_ATTENUATION_AB, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy,&
epsilondev_xz,epsilondev_yz,epsilon_trace_over_3, &
ANISOTROPY,NSPEC_ANISO, &
- c11store,c12store,c13store,c14store,c15store,c16store,&
- c22store,c23store,c24store,c25store,c26store,c33store,&
- c34store,c35store,c36store,c44store,c45store,c46store,&
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
c55store,c56store,c66store, &
SIMULATION_TYPE,COMPUTE_AND_STORE_STRAIN,NSPEC_STRAIN_ONLY, &
- NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ADJOINT,&
+ NSPEC_BOUN,NSPEC2D_MOHO,NSPEC_ADJOINT, &
is_moho_top,is_moho_bot, &
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
- num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic,&
+ num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
phase_ispec_inner_elastic )
use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM, &
@@ -61,9 +61,12 @@
integer :: NSPEC_AB,NGLOB_AB
-! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+! displacement, velocity and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
+! time step
+ real(kind=CUSTOM_REAL) :: deltat
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -145,6 +148,14 @@
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL) tempx1l_att,tempx2l_att,tempx3l_att
+ real(kind=CUSTOM_REAL) tempy1l_att,tempy2l_att,tempy3l_att
+ real(kind=CUSTOM_REAL) tempz1l_att,tempz2l_att,tempz3l_att
+
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att;
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att;
+
! local anisotropy parameters
real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
@@ -219,6 +230,47 @@
tempz3l = tempz3l + displ(3,iglob)*hp3
enddo
+ if( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
+ tempx1l_att = tempx1l
+ tempx2l_att = tempx2l
+ tempx3l_att = tempx3l
+
+ tempy1l_att = tempy1l
+ tempy2l_att = tempy2l
+ tempy3l_att = tempy3l
+
+ tempz1l_att = tempz1l
+ tempz2l_att = tempz2l
+ tempz3l_att = tempz3l
+
+ ! use first order Taylor expansion of displacement for local storage of stresses
+ ! at this current time step, to fix attenuation in a consistent way
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1l_att = tempx1l_att + deltat*veloc(1,iglob)*hp1
+ tempy1l_att = tempy1l_att + deltat*veloc(2,iglob)*hp1
+ tempz1l_att = tempz1l_att + deltat*veloc(3,iglob)*hp1
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2l_att = tempx2l_att + deltat*veloc(1,iglob)*hp2
+ tempy2l_att = tempy2l_att + deltat*veloc(2,iglob)*hp2
+ tempz2l_att = tempz2l_att + deltat*veloc(3,iglob)*hp2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3l_att = tempx3l_att + deltat*veloc(1,iglob)*hp3
+ tempy3l_att = tempy3l_att + deltat*veloc(2,iglob)*hp3
+ tempz3l_att = tempz3l_att + deltat*veloc(3,iglob)*hp3
+ enddo
+ endif
+
! get derivatives of ux, uy and uz with respect to x, y and z
xixl = xix(i,j,k,ispec)
xiyl = xiy(i,j,k,ispec)
@@ -276,15 +328,43 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! computes deviatoric strain attenuation and/or for kernel calculations
- if (COMPUTE_AND_STORE_STRAIN) then
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
- epsilondev_xx_loc(i,j,k) = duxdxl - templ
- epsilondev_yy_loc(i,j,k) = duydyl - templ
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att
+ duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att
+ duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att
+
+ duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att
+ duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att
+ duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att
+
+ duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att
+ duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att
+ duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att
+
+ ! precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
+
+ ! compute deviatoric strain
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ epsilondev_xx_loc(i,j,k) = duxdxl_att - epsilon_trace_over_3(i,j,k,ispec)
+ epsilondev_yy_loc(i,j,k) = duydyl_att - epsilon_trace_over_3(i,j,k,ispec)
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! computes deviatoric strain attenuation and/or for kernel calculations
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( SIMULATION_TYPE == 3 ) epsilon_trace_over_3(i,j,k,ispec) = templ
+ epsilondev_xx_loc(i,j,k) = duxdxl - templ
+ epsilondev_yy_loc(i,j,k) = duydyl - templ
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
kappal = kappastore(i,j,k,ispec)
More information about the CIG-COMMITS
mailing list