[cig-commits] r8429 - seismo/2D/SPECFEM2D/trunk/SPECFEM90
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:45:34 PST 2007
Author: walter
Date: 2007-12-07 15:45:33 -0800 (Fri, 07 Dec 2007)
New Revision: 8429
Added:
seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_attenuation.f90
Modified:
seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile
seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h
seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
Log:
finished implementation of attenuation in 2D code
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile 2004-09-11 10:14:06 UTC (rev 8428)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile 2007-12-07 23:45:33 UTC (rev 8429)
@@ -25,7 +25,7 @@
EXEC = xspecfem2D
OBJS = $O/checkgrid.o $O/datim.o $O/defarrays.o\
$O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o $O/define_derivative_matrices.o\
- $O/plotpost.o $O/positrec.o $O/positsource.o $O/q49spec.o\
+ $O/plotpost.o $O/positrec.o $O/positsource.o $O/q49spec.o $O/compute_gradient_attenuation.o\
$O/specfem2D.o $O/write_seismograms.o $O/createnum_fast.o $O/createnum_slow.o $O/q49shape.o
default: all
@@ -84,6 +84,9 @@
$O/specfem2D.o: specfem2D.f90 constants.h
${F90} $(FLAGS) -c -o $O/specfem2D.o specfem2D.f90
+$O/compute_gradient_attenuation.o: compute_gradient_attenuation.f90 constants.h
+ ${F90} $(FLAGS) -c -o $O/compute_gradient_attenuation.o compute_gradient_attenuation.f90
+
$O/write_seismograms.o: write_seismograms.f90 constants.h
${F90} $(FLAGS) -c -o $O/write_seismograms.o write_seismograms.f90
Added: seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_attenuation.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_attenuation.f90 2004-09-11 10:14:06 UTC (rev 8428)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_attenuation.f90 2007-12-07 23:45:33 UTC (rev 8429)
@@ -0,0 +1,88 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.0
+! ------------------------------
+!
+! Dimitri Komatitsch
+! Universite de Pau et des Pays de l'Adour, France
+!
+! (c) May 2004
+!
+!========================================================================
+
+ subroutine compute_gradient_attenuation(displ,duxdxl,duzdxl,duxdzl,duzdzl,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
+
+! compute Grad(displ) for attenuation
+
+ implicit none
+
+ include "constants.h"
+
+ integer NSPEC,npoin
+
+ integer, dimension(NGLLX,NGLLZ,NSPEC) :: ibool
+
+ double precision, dimension(NGLLX,NGLLZ,NSPEC) :: duxdxl,duzdxl,duxdzl,duzdzl,xix,xiz,gammax,gammaz
+
+ double precision, dimension(NDIME,npoin) :: displ
+
+! array with derivatives of Lagrange polynomials
+ double precision, dimension(NGLLX,NGLLX) :: hprime_xx
+ double precision, dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+! local variables
+ integer i,j,k,ispec,iglob
+
+! space derivatives
+ double precision tempx1l,tempx2l,tempz1l,tempz2l
+ double precision hp1,hp2
+
+! jacobian
+ double precision xixl,xizl,gammaxl,gammazl
+
+! loop over spectral elements
+ do ispec = 1,NSPEC
+
+! double loop over GLL to compute and store gradients
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+! derivative along x
+ tempx1l = ZERO
+ tempz1l = ZERO
+ do k = 1,NGLLX
+ hp1 = hprime_xx(k,i)
+ iglob = ibool(k,j,ispec)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempz1l = tempz1l + displ(2,iglob)*hp1
+ enddo
+
+! derivative along z
+ tempx2l = ZERO
+ tempz2l = ZERO
+ do k = 1,NGLLZ
+ hp2 = hprime_zz(k,j)
+ iglob = ibool(i,k,ispec)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempz2l = tempz2l + displ(2,iglob)*hp2
+ enddo
+
+ xixl = xix(i,j,ispec)
+ xizl = xiz(i,j,ispec)
+ gammaxl = gammax(i,j,ispec)
+ gammazl = gammaz(i,j,ispec)
+
+! derivatives of displacement
+ duxdxl(i,j,ispec) = tempx1l*xixl + tempx2l*gammaxl
+ duxdzl(i,j,ispec) = tempx1l*xizl + tempx2l*gammazl
+
+ duzdxl(i,j,ispec) = tempz1l*xixl + tempz2l*gammaxl
+ duzdzl(i,j,ispec) = tempz1l*xizl + tempz2l*gammazl
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine compute_gradient_attenuation
+
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h 2004-09-11 10:14:06 UTC (rev 8428)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h 2007-12-07 23:45:33 UTC (rev 8429)
@@ -87,9 +87,6 @@
! anisotropic copper crystal (cubic symmetry)
!
-! switch anisotropy on or off
- logical, parameter :: TURN_ANISOTROPY_ON = .false.
-
! regular c_ijkl with no rotation
double precision, parameter :: c11val = 169.d9
double precision, parameter :: c12val = 122.d9
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2004-09-11 10:14:06 UTC (rev 8428)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2007-12-07 23:45:33 UTC (rev 8429)
@@ -20,7 +20,7 @@
! version 5.0 : - got rid of useless routines, suppressed commons etc.
! - weak formulation based explicitly on stress tensor
! - implementation of full anisotropy
-! - implementation of attenuation
+! - implementation of attenuation based on memory variables
! version 5.0 is based on SPECFEM2D version 4.2
! (c) June 1998 by Dimitri Komatitsch, Harvard University, USA
@@ -80,7 +80,8 @@
double precision xixl,xizl,gammaxl,gammazl,jacobianl
! material properties of the elastic medium
- double precision mul,lambdal,lambdalplus2mul
+ double precision mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed
+ double precision mul_unrelaxed,lambdal_unrelaxed,lambdalplus2mul_unrelaxed
double precision, dimension(:), allocatable :: xirec,etarec
@@ -107,13 +108,13 @@
double precision f0,t0,factor,a,angleforce,ricker,displnorm_all
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
- lambdalSmin,lambdalSmax,lambdalPmin,lambdalPmax,vpmin,vpmax
+ lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax,vpmin,vpmax
integer icolor,inumber,isubsamp,ivecttype,itaff,nrec,isismostype
integer numat,ngnod,nspec,iptsdisp,nelemabs
logical interpol,imeshvect,imodelvect,iboundvect,ireadmodel,initialfield, &
- ioutputgrid,ignuplot
+ ioutputgrid,ignuplot,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
double precision cutvect,anglerec
@@ -124,10 +125,12 @@
logical, dimension(:,:), allocatable :: codeabs
-! DK DK feb99 for attenuation
+! for attenuation
+ integer nspec_allocate
+
double precision tau_epsilon_nu1_mech1,tau_sigma_nu1_mech1, &
- tau_epsilon_nu2_mech1,tau_sigma_nu2_mech1,tau_epsilon_nu1_mech2, &
- tau_sigma_nu1_mech2,tau_epsilon_nu2_mech2,tau_sigma_nu2_mech2
+ tau_epsilon_nu2_mech1,tau_sigma_nu2_mech1,tau_epsilon_nu1_mech2, &
+ tau_sigma_nu1_mech2,tau_epsilon_nu2_mech2,tau_sigma_nu2_mech2
double precision Un,Unp1,tauinv,Sn,Snp1,theta_n,theta_np1
double precision phi_nu1_mech1,phi_nu2_mech1,phi_nu1_mech2,phi_nu2_mech2
@@ -141,7 +144,7 @@
double precision, dimension(:,:,:), allocatable :: &
e1_mech1,e11_mech1,e13_mech1,e1_mech2,e11_mech2,e13_mech2, &
- DxUx_n,DzUz_n,DxUz_n,DzUx_n,DxUx_np1,DzUz_np1,DxUz_np1,DzUx_np1
+ duxdxl_n,duzdzl_n,duzdxl_n,duxdzl_n,duxdxl_np1,duzdzl_np1,duzdxl_np1,duxdzl_np1
! title of the plot
character(len=60) stitle
@@ -207,13 +210,13 @@
read(IIN,*) isismostype,ivecttype
read(IIN,40) datlin
- read(IIN,*) ireadmodel,ioutputgrid
+ read(IIN,*) ireadmodel,ioutputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
!---- check parameters read
write(IOUT,200) npgeo,NDIME
write(IOUT,600) itaff,icolor,inumber
write(IOUT,700) nrec,isismostype,anglerec
- write(IOUT,750) initialfield,ireadmodel,ioutputgrid
+ write(IOUT,750) initialfield,ireadmodel,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,ioutputgrid
write(IOUT,800) ivecttype,100.d0*cutvect,isubsamp
!---- read time step
@@ -317,29 +320,30 @@
allocate(knods(ngnod,nspec))
allocate(ibool(NGLLX,NGLLZ,nspec))
-! DK DK feb99 for attenuation
-!! DK DK declarer le flag TURN_ATTGEN_ON et le lire depuis le parfile
-!! DK DK aussi faire attention a ne pas mettre aniso et attenuation (tester et stop si les deux)
-!! DK DK aussi mettre a jour comments feb99 -> 2004
-if(TURN_ATTGEN_ON) then
- allocate(e1_mech1(NGLLX,NGLLZ,nspec))
- allocate(e11_mech1(NGLLX,NGLLZ,nspec))
- allocate(e13_mech1(NGLLX,NGLLZ,nspec))
- allocate(e1_mech2(NGLLX,NGLLZ,nspec))
- allocate(e11_mech2(NGLLX,NGLLZ,nspec))
- allocate(e13_mech2(NGLLX,NGLLZ,nspec))
- allocate(DxUx_n(NGLLX,NGLLZ,nspec))
- allocate(DzUz_n(NGLLX,NGLLZ,nspec))
- allocate(DxUz_n(NGLLX,NGLLZ,nspec))
- allocate(DzUx_n(NGLLX,NGLLZ,nspec))
- allocate(DxUx_np1(NGLLX,NGLLZ,nspec))
- allocate(DzUz_np1(NGLLX,NGLLZ,nspec))
- allocate(DxUz_np1(NGLLX,NGLLZ,nspec))
- allocate(DzUx_np1(NGLLX,NGLLZ,nspec))
-else
- allocate_size_one
-endif
+! for attenuation
+ if(TURN_ANISOTROPY_ON .and. TURN_ATTENUATION_ON) stop 'cannot have anisotropy and attenuation both turned on in current version'
+ if(TURN_ATTENUATION_ON) then
+ nspec_allocate = nspec
+ else
+ nspec_allocate = 1
+ endif
+
+ allocate(e1_mech1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(e11_mech1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(e13_mech1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(e1_mech2(NGLLX,NGLLZ,nspec_allocate))
+ allocate(e11_mech2(NGLLX,NGLLZ,nspec_allocate))
+ allocate(e13_mech2(NGLLX,NGLLZ,nspec_allocate))
+ allocate(duxdxl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(duzdzl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(duzdxl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(duxdzl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(duxdxl_np1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(duzdzl_np1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(duzdxl_np1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(duxdzl_np1(NGLLX,NGLLZ,nspec_allocate))
+
! --- allocate arrays for absorbing boundary conditions
if(nelemabs <= 0) then
nelemabs = 1
@@ -526,8 +530,8 @@
call defarrays(vpext,vsext,rhoext,density,elastcoef, &
xigll,yigll,xix,xiz,gammax,gammaz,a11,a12, &
ibool,kmato,coord,gltfu,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,lambdalSmin,lambdalSmax, &
- lambdalPmin,lambdalPmax,vpmin,vpmax,ireadmodel,nspec,numat)
+ cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax, &
+ lambdal_Pmin,lambdal_Pmax,vpmin,vpmax,ireadmodel,nspec,numat)
! build the global mass matrix once and for all
rmass(:) = ZERO
@@ -553,7 +557,7 @@
!---- verifier le maillage, la stabilite et le nb de points par lambda
!
call checkgrid(deltat,gltfu,initialfield,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,lambdalSmin,lambdalSmax,lambdalPmin,lambdalPmax)
+ cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax)
!
!---- initialiser sismogrammes
@@ -597,10 +601,11 @@
print *,'Max norm of initial displacement = ',maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
endif
-! DK DK attenuation constants from Carcione 1993 Geophysics volume 58 pages 111 and 112
+! attenuation constants from Carcione 1993 Geophysics volume 58 pages 111 and 112
! for two memory-variables mechanisms.
! beware: these values implement specific values of the quality factor Q,
! see Carcione 1993 for details
+
tau_epsilon_nu1_mech1 = 0.0334d0
tau_sigma_nu1_mech1 = 0.0303d0
tau_epsilon_nu2_mech1 = 0.0352d0
@@ -610,20 +615,18 @@
tau_epsilon_nu2_mech2 = 0.0029d0
tau_sigma_nu2_mech2 = 0.0024d0
- inv_tau_sigma_nu1_mech1 = 1.d0 / tau_sigma_nu1_mech1
- inv_tau_sigma_nu2_mech1 = 1.d0 / tau_sigma_nu2_mech1
- inv_tau_sigma_nu1_mech2 = 1.d0 / tau_sigma_nu1_mech2
- inv_tau_sigma_nu2_mech2 = 1.d0 / tau_sigma_nu2_mech2
+ inv_tau_sigma_nu1_mech1 = ONE / tau_sigma_nu1_mech1
+ inv_tau_sigma_nu2_mech1 = ONE / tau_sigma_nu2_mech1
+ inv_tau_sigma_nu1_mech2 = ONE / tau_sigma_nu1_mech2
+ inv_tau_sigma_nu2_mech2 = ONE / tau_sigma_nu2_mech2
- phi_nu1_mech1 = (1.d0 - tau_epsilon_nu1_mech1/tau_sigma_nu1_mech1) / tau_sigma_nu1_mech1
- phi_nu2_mech1 = (1.d0 - tau_epsilon_nu2_mech1/tau_sigma_nu2_mech1) / tau_sigma_nu2_mech1
- phi_nu1_mech2 = (1.d0 - tau_epsilon_nu1_mech2/tau_sigma_nu1_mech2) / tau_sigma_nu1_mech2
- phi_nu2_mech2 = (1.d0 - tau_epsilon_nu2_mech2/tau_sigma_nu2_mech2) / tau_sigma_nu2_mech2
+ phi_nu1_mech1 = (ONE - tau_epsilon_nu1_mech1/tau_sigma_nu1_mech1) / tau_sigma_nu1_mech1
+ phi_nu2_mech1 = (ONE - tau_epsilon_nu2_mech1/tau_sigma_nu2_mech1) / tau_sigma_nu2_mech1
+ phi_nu1_mech2 = (ONE - tau_epsilon_nu1_mech2/tau_sigma_nu1_mech2) / tau_sigma_nu1_mech2
+ phi_nu2_mech2 = (ONE - tau_epsilon_nu2_mech2/tau_sigma_nu2_mech2) / tau_sigma_nu2_mech2
- Mu_nu1 = 1.d0 - (1.d0 - tau_epsilon_nu1_mech1/tau_sigma_nu1_mech1) &
- - (1.d0 - tau_epsilon_nu1_mech2/tau_sigma_nu1_mech2)
- Mu_nu2 = 1.d0 - (1.d0 - tau_epsilon_nu2_mech1/tau_sigma_nu2_mech1) &
- - (1.d0 - tau_epsilon_nu2_mech2/tau_sigma_nu2_mech2)
+ Mu_nu1 = ONE - (ONE - tau_epsilon_nu1_mech1/tau_sigma_nu1_mech1) - (ONE - tau_epsilon_nu1_mech2/tau_sigma_nu1_mech2)
+ Mu_nu2 = ONE - (ONE - tau_epsilon_nu2_mech1/tau_sigma_nu2_mech1) - (ONE - tau_epsilon_nu2_mech2/tau_sigma_nu2_mech2)
deltatsquare = deltat * deltat
deltatcube = deltatsquare * deltat
@@ -653,8 +656,9 @@
endif
endif
-! DK DK calculer Grad(U) a l'instant n pour attenuation
- call calc_gradient(DxUx_n,DzUz_n,DxUz_n,DzUx_n,xjaci,hprime,hTprime,ibool,displ,NGLLX,npoin,ndime,nspec)
+! compute Grad(displ) at time step n for attenuation
+ if(TURN_ATTENUATION_ON) call compute_gradient_attenuation(displ,duxdxl_n,duzdxl_n, &
+ duxdzl_n,duzdzl_n,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
! `predictor' update displacement using finite-difference time scheme (Newmark)
displ(:,:) = displ(:,:) + deltat*veloc(:,:) + deltatsquareover2*accel(:,:)
@@ -664,10 +668,10 @@
! integration over spectral elements
do ispec = 1,NSPEC
-! get elastic parameters of current spectral element
- lambdal = elastcoef(1,kmato(ispec))
- mul = elastcoef(2,kmato(ispec))
- lambdalplus2mul = lambdal + 2.d0*mul
+! get relaxed elastic parameters of current spectral element
+ lambdal_relaxed = elastcoef(1,kmato(ispec))
+ mul_relaxed = elastcoef(2,kmato(ispec))
+ lambdalplus2mul_relaxed = lambdal_relaxed + TWO*mul_relaxed
! first double loop over GLL to compute and store gradients
do j = 1,NGLLZ
@@ -679,11 +683,16 @@
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
- mul = rhol*csl*csl
- lambdal = rhol*cpl*cpl - 2.d0*mul
- lambdalplus2mul = lambdal + 2.d0*mul
+ mul_relaxed = rhol*csl*csl
+ lambdal_relaxed = rhol*cpl*cpl - TWO*mul_relaxed
+ lambdalplus2mul_relaxed = lambdal_relaxed + TWO*mul_relaxed
endif
+! compute unrelaxed elastic coefficients from formulas in Carcione 1993 page 111
+ lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1 - mul_relaxed * Mu_nu2
+ mul_unrelaxed = mul_relaxed * Mu_nu2
+ lambdalplus2mul_unrelaxed = lambdal_unrelaxed + TWO*mul_unrelaxed
+
! derivative along x
tempx1l = ZERO
tempz1l = ZERO
@@ -716,71 +725,30 @@
duzdxl = tempz1l*xixl + tempz2l*gammaxl
duzdzl = tempz1l*xizl + tempz2l*gammazl
+! compute stress tensor (include attenuation or anisotropy if needed)
-
-
-
-
-
-
-
-
-
-! compute stress tensor
-
-! DK DK for attenuation
-
if(TURN_ATTENUATION_ON) then
! compute the stress using the unrelaxed Lame parameters (Carcione page 111)
- sigma_xx = (rlambda_unrelax + 2.d0*rmu_unrelax)*DxUx + rlambda_unrelax*DzUz
- sigma_zz = rlambda_unrelax*DxUx + (rlambda_unrelax + 2.d0*rmu_unrelax)*DzUz
- sigma_xz = rmu_unrelax*(DxUz + DzUx)
+ sigma_xx = lambdalplus2mul_unrelaxed*duxdxl + lambdal_unrelaxed*duzdzl
+ sigma_xz = mul_unrelaxed*(duzdxl + duxdzl)
+ sigma_zz = lambdalplus2mul_unrelaxed*duzdzl + lambdal_unrelaxed*duxdxl
! add the memory variables using the relaxed parameters (Carcione page 111)
! beware: there is a bug in Carcione's equation for sigma_zz
- sigma_xx = sigma_xx + (rlambda_relax + rmu_relax)* (e1_mech1(i,j,k) + e1_mech2(i,j,k)) &
- + 2.d0 * rmu_relax * (e11_mech1(i,j,k) + e11_mech2(i,j,k))
- sigma_zz = sigma_zz + (rlambda_relax + rmu_relax)* (e1_mech1(i,j,k) + e1_mech2(i,j,k)) &
- - 2.d0 * rmu_relax * (e11_mech1(i,j,k) + e11_mech2(i,j,k))
- sigma_xz = sigma_xz + rmu_relax * (e13_mech1(i,j,k) + e13_mech2(i,j,k))
+ sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed)* (e1_mech1(i,j,k) + e1_mech2(i,j,k)) + TWO * mul_relaxed * (e11_mech1(i,j,k) + e11_mech2(i,j,k))
+ sigma_xz = sigma_xz + mul_relaxed * (e13_mech1(i,j,k) + e13_mech2(i,j,k))
+ sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed)* (e1_mech1(i,j,k) + e1_mech2(i,j,k)) - TWO * mul_relaxed * (e11_mech1(i,j,k) + e11_mech2(i,j,k))
else
-! DK DK no attenuation
+! no attenuation
+ sigma_xx = lambdalplus2mul_relaxed*duxdxl + lambdal_relaxed*duzdzl
+ sigma_xz = mul_relaxed*(duzdxl + duxdzl)
+ sigma_zz = lambdalplus2mul_relaxed*duzdzl + lambdal_relaxed*duxdxl
-! old clean
- sigma_xx = lambdalplus2mul*duxdxl + lambdal*duzdzl
- sigma_xz = mul*(duzdxl + duxdzl)
- sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl
-
-! new dirty
- sigma_xx = (rlambda_relax + 2.d0*rmu_relax)*DxUx + rlambda_relax*DzUz
- sigma_zz = rlambda_relax*DxUx + (rlambda_relax + 2.d0*rmu_relax)*DzUz
- sigma_xz = rmu_relax*(DxUz + DzUx)
-
endif
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
! full anisotropy
if(TURN_ANISOTROPY_ON) then
@@ -880,12 +848,12 @@
ispec = numabs(ispecabs)
! get elastic parameters of current spectral element
- lambdal = elastcoef(1,kmato(ispec))
- mul = elastcoef(2,kmato(ispec))
+ lambdal_relaxed = elastcoef(1,kmato(ispec))
+ mul_relaxed = elastcoef(2,kmato(ispec))
rhol = density(kmato(ispec))
- rKvol = lambdal + 2.d0*mul/3.d0
- cpl = dsqrt((rKvol + 4.d0*mul/3.d0)/rhol)
- csl = dsqrt(mul/rhol)
+ rKvol = lambdal_relaxed + TWO*mul_relaxed/3.d0
+ cpl = dsqrt((rKvol + 4.d0*mul_relaxed/3.d0)/rhol)
+ csl = dsqrt(mul_relaxed/rhol)
!--- left absorbing boundary
@@ -909,8 +877,8 @@
rho_vp = rhol*cpl
rho_vs = rhol*csl
- nx = -1.d0
- nz = 0.d0
+ nx = -ONE
+ nz = ZERO
vx = veloc(1,iglob)
vz = veloc(2,iglob)
@@ -950,8 +918,8 @@
rho_vp = rhol*cpl
rho_vs = rhol*csl
- nx = 1.d0
- nz = 0.d0
+ nx = ONE
+ nz = ZERO
vx = veloc(1,iglob)
vz = veloc(2,iglob)
@@ -997,8 +965,8 @@
rho_vp = rhol*cpl
rho_vs = rhol*csl
- nx = 0.d0
- nz = -1.d0
+ nx = ZERO
+ nz = -ONE
vx = veloc(1,iglob)
vz = veloc(2,iglob)
@@ -1044,8 +1012,8 @@
rho_vp = rhol*cpl
rho_vs = rhol*csl
- nx = 0.d0
- nz = 1.d0
+ nx = ZERO
+ nz = ONE
vx = veloc(1,iglob)
vz = veloc(2,iglob)
@@ -1079,7 +1047,7 @@
! Ricker wavelet for the source time function
a = pi*pi*f0*f0
- ricker = - factor * (1.d0-2.d0*a*(time-t0)**2)*exp(-a*(time-t0)**2)
+ ricker = - factor * (ONE-TWO*a*(time-t0)**2)*exp(-a*(time-t0)**2)
! --- collocated force
if(nint(gltfu(2)) == 1) then
@@ -1111,25 +1079,20 @@
! `corrector' update velocity
veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
+! implement attenuation
+ if(TURN_ATTENUATION_ON) then
+! compute Grad(displ) at time step n+1 for attenuation
+ call compute_gradient_attenuation(displ,duxdxl_np1,duzdxl_np1, &
+ duxdzl_np1,duzdzl_np1,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
-
-
-
-
-
-! DK DK feb99 mettre a jour memory variables par Runge-Kutta ordre 4 ici for attenuation
-
-
-! calculer Grad(U) a l'instant n+1 pour attenuation
- call calc_gradient(DxUx_np1,DzUz_np1,DxUz_np1,DzUx_np1,xjaci,hprime,hTprime,ibool,displ,NGLLX,npoin,ndime,nspec)
-
+! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
do k=1,nspec
do j=1,NGLLZ
do i=1,NGLLX
- theta_n = DxUx_n(i,j,k) + DzUz_n(i,j,k)
- theta_np1 = DxUx_np1(i,j,k) + DzUz_np1(i,j,k)
+ theta_n = duxdxl_n(i,j,k) + duzdzl_n(i,j,k)
+ theta_np1 = duxdxl_np1(i,j,k) + duzdzl_np1(i,j,k)
! evolution e1_mech1
Un = e1_mech1(i,j,k)
@@ -1165,8 +1128,8 @@
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = (DxUx_n(i,j,k) - theta_n/2.d0) * phi_nu2_mech1
- Snp1 = (DxUx_np1(i,j,k) - theta_np1/2.d0) * phi_nu2_mech1
+ Sn = (duxdxl_n(i,j,k) - theta_n/TWO) * phi_nu2_mech1
+ Snp1 = (duxdxl_np1(i,j,k) - theta_np1/TWO) * phi_nu2_mech1
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
@@ -1179,8 +1142,8 @@
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = (DxUx_n(i,j,k) - theta_n/2.d0) * phi_nu2_mech2
- Snp1 = (DxUx_np1(i,j,k) - theta_np1/2.d0) * phi_nu2_mech2
+ Sn = (duxdxl_n(i,j,k) - theta_n/TWO) * phi_nu2_mech2
+ Snp1 = (duxdxl_np1(i,j,k) - theta_np1/TWO) * phi_nu2_mech2
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
@@ -1193,8 +1156,8 @@
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = (DzUx_n(i,j,k) + DxUz_n(i,j,k)) * phi_nu2_mech1
- Snp1 = (DzUx_np1(i,j,k) + DxUz_np1(i,j,k)) * phi_nu2_mech1
+ Sn = (duxdzl_n(i,j,k) + duzdxl_n(i,j,k)) * phi_nu2_mech1
+ Snp1 = (duxdzl_np1(i,j,k) + duzdxl_np1(i,j,k)) * phi_nu2_mech1
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
@@ -1207,8 +1170,8 @@
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = (DzUx_n(i,j,k) + DxUz_n(i,j,k)) * phi_nu2_mech2
- Snp1 = (DzUx_np1(i,j,k) + DxUz_np1(i,j,k)) * phi_nu2_mech2
+ Sn = (duxdzl_n(i,j,k) + duzdxl_n(i,j,k)) * phi_nu2_mech2
+ Snp1 = (duxdzl_np1(i,j,k) + duzdxl_np1(i,j,k)) * phi_nu2_mech2
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
@@ -1219,22 +1182,8 @@
enddo
enddo
-! DK DK feb99 mettre a jour memory variables par Runge-Kutta ordre 4 ici for attenuation
+ endif ! end of test on attenuation
-
-
-
-
-
-
-
-
-
-
-
-
-
-
!---- display max of norm of displacement
if(mod(it,itaff) == 0) then
displnorm_all = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
@@ -1359,6 +1308,8 @@
750 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
'Read external initial field or not . .(initialfield) = ',l6/5x, &
'Read external velocity model or not. . .(ireadmodel) = ',l6/5x, &
+ 'Turn anisotropy on or off. . . .(TURN_ANISOTROPY_ON) = ',l6/5x, &
+ 'Turn attenuation on or off. . .(TURN_ATTENUATION_ON) = ',l6/5x, &
'Save grid in external file or not . . .(ioutputgrid) = ',l6)
800 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
'Vector display type . . . . . . . . . . .(ivecttype) = ',i6/5x, &
More information about the cig-commits
mailing list