[cig-commits] r22906 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon Sep 30 18:02:11 PDT 2013
Author: xie.zhinan
Date: 2013-09-30 18:02:10 -0700 (Mon, 30 Sep 2013)
New Revision: 22906
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
remove ADE formualtion of viscoelastic wave equations using Newmark scheme to update memory variables. Add convolutional formulation of viscoelastic formulation and using second order convolutional scheme to update memory variables. Clean the code asscociated with viscoelastic simulation
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90 2013-09-30 17:22:41 UTC (rev 22905)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90 2013-10-01 01:02:10 UTC (rev 22906)
@@ -47,15 +47,12 @@
initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
accelw_poroelastic,velocw_poroelastic,displw_poroelastic,velocs_poroelastic,displs_poroelastic,&
- b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
+ displs_poroelastic_old,b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
- jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
- e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu2,&
- phi_nu2,Mu_nu2,N_SLS, &
- rx_viscous,rz_viscous,theta_e,theta_s,&
- b_viscodampx,b_viscodampz,&
+ jacobian,source_time_function,sourcearray,adj_sourcearrays, &
+ e11,e13,hprime_xx,hprimewgll_xx,hprime_zz,hprimewgll_zz,wxgll,wzgll,&
+ inv_tau_sigma_nu2,phi_nu2,Mu_nu2,N_SLS, &
+ rx_viscous,rz_viscous,theta_e,theta_s,b_viscodampx,b_viscodampz,&
ibegin_edge1_poro,iend_edge1_poro,ibegin_edge3_poro,iend_edge3_poro, &
ibegin_edge4_poro,iend_edge4_poro,ibegin_edge2_poro,iend_edge2_poro,&
C_k,M_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
@@ -96,7 +93,7 @@
logical, dimension(4,nelemabs) :: codeabs
real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: accelw_poroelastic,velocw_poroelastic,displw_poroelastic,&
- displs_poroelastic,velocs_poroelastic
+ displs_poroelastic,displs_poroelastic_old,velocs_poroelastic
real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic
double precision, dimension(2,numat) :: density
double precision, dimension(3,numat) :: permeability
@@ -115,8 +112,6 @@
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11,e13
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_veloc,e13_veloc
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_accel,e13_accel
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_LDDRK,e13_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS,stage_time_scheme) :: e11_force_RK, e13_force_RK
@@ -125,8 +120,9 @@
real(kind=CUSTOM_REAL) :: e11_sum,e13_sum
integer :: i_sls
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
- dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) ::dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ !nsub1 denote discrete time step n-1
+ dux_dxl_nsub1,duz_dzl_nsub1,duz_dxl_nsub1,dux_dzl_nsub1
! viscous attenuation
double precision, dimension(NGLLX,NGLLZ,nspec) :: rx_viscous
@@ -143,7 +139,6 @@
real(kind=CUSTOM_REAL), dimension(NGLLX) :: wxgll
real(kind=CUSTOM_REAL), dimension(NGLLZ) :: wzgll
-!
double precision :: f0,freq0,Q0,w_c
! Parameter for LDDRK time scheme
@@ -195,132 +190,107 @@
real(kind=CUSTOM_REAL) :: cpIsquare,cpIIsquare,cssquare,cpIl,cpIIl,csl
! for attenuation
- real(kind=CUSTOM_REAL) :: phinu2,tauinvnu2,theta_n_u,theta_n_v
+ real(kind=CUSTOM_REAL) :: phinu2,tauinvnu2,theta_n_u,theta_nsub1_u
+ real(kind=CUSTOM_REAL) :: bb,coef0,coef1,coef2
! implement attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
- ! compute Grad(displs_poroelastic) at time step n for attenuation
+! compute Grad(displs_poroelastic) at time step n for attenuation
call compute_gradient_attenuation(displs_poroelastic,dux_dxl_n,duz_dxl_n, &
- dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,nglob)
+ dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,nglob)
- ! compute Grad(velocs_poroelastic) at time step n for attenuation
- call compute_gradient_attenuation(velocs_poroelastic,dvx_dxl_n,dvz_dxl_n, &
- dvx_dzl_n,dvz_dzl_n,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,nglob)
+! compute Grad(displs_poroelastic) at time step n-1 for attenuation
+ call compute_gradient_attenuation(displs_poroelastic_old,dux_dxl_nsub1,duz_dxl_nsub1, &
+ dux_dzl_nsub1,duz_dzl_nsub1,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,nglob)
! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
! loop over spectral elements
do ispec = 1,nspec
if (poroelastic(ispec)) then
-
- do j=1,NGLLZ
- do i=1,NGLLX
+ do j=1,NGLLZ; do i=1,NGLLX
theta_n_u = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
- theta_n_v = dvx_dxl_n(i,j,ispec) + dvz_dzl_n(i,j,ispec)
+ theta_nsub1_u = dux_dxl_nsub1(i,j,ispec) + duz_dzl_nsub1(i,j,ispec)
-! loop on all the standard linear solids
+ ! loop on all the standard linear solids
do i_sls = 1,N_SLS
+ phinu2 = phi_nu2(i,j,ispec,i_sls)
+ tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
+ ! update e1, e11, e13 in convolution formation with modified recursive convolution scheme on basis of
+ ! second-order accurate convolution term calculation from equation (21) of
+ ! Shumin Wang, Robert Lee, and Fernando L. Teixeira,
+ ! Anisotropic-Medium PML for Vector FETD With Modified Basis Functions,
+ ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006)
+ ! evolution e1 ! no need since we are just considering shear attenuation
+ if(stage_time_scheme == 1) then
+ bb = tauinvnu2; coef0 = exp(-bb * deltat)
+ if ( abs(bb) > 1e-5_CUSTOM_REAL ) then
+ coef1 = (1._CUSTOM_REAL - exp(-bb * deltat / 2._CUSTOM_REAL)) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-bb* deltat / 2._CUSTOM_REAL)) * exp(-bb * deltat / 2._CUSTOM_REAL)/ bb
+ else
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
+ endif
-! evolution e1 ! no need since we are just considering shear attenuation
-! if(stage_time_scheme == 1) then
-! Un = e11(i,j,ispec,i_sls)
-! tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
-! tauinvsquare = tauinv * tauinv
-! tauinvcube = tauinvsquare * tauinv
-! tauinvUn = tauinv * Un
-! Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i,j,ispec,i_sls)
-! Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i,j,ispec,i_sls)
-! Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
-! twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
-! fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
-! deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
-! e11(i,j,ispec,i_sls) = Unp1
-! endif
+ e11(i,j,ispec,i_sls) = coef0 * e11(i,j,ispec,i_sls) + &
+ phinu2 * (coef1 * (dux_dxl_n(i,j,ispec)-theta_n_u/TWO) + &
+ coef2 * (dux_dxl_nsub1(i,j,ispec)-theta_nsub1_u/TWO))
- if(stage_time_scheme == 1) then
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) + deltat*e11_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e11_accel(i,j,ispec,i_sls)
- e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- e11_accel(i,j,ispec,i_sls) = ((dvx_dxl_n(i,j,ispec)-theta_n_v/TWO) * phinu2- &
- e11_veloc(i,j,ispec,i_sls)*tauinvnu2) &
- /(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
- e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- endif
+ e13(i,j,ispec,i_sls) = coef0 * e13(i,j,ispec,i_sls) + &
+ phinu2 * (coef1 * (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) + &
+ coef2 * (dux_dzl_nsub1(i,j,ispec) + duz_dxl_nsub1(i,j,ispec)))
+ endif
+ ! update e1, e11, e13 in ADE formation with LDDRK scheme
+ ! evolution e1 ! no need since we are just considering shear attenuation
+ if(stage_time_scheme == 6) then
+ e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) + &
+ deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) - &
+ deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
+ e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*e11_LDDRK(i,j,ispec,i_sls)
- if(stage_time_scheme == 6) then
- e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
- + deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) &
- - deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*e11_LDDRK(i,j,ispec,i_sls)
- endif
+ e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) + &
+ deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) - &
+ deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
+ e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)
+ endif
- if(stage_time_scheme == 4) then
- e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2- &
- e11(i,j,ispec,i_sls) * tauinvnu2)
+ ! update e1, e11, e13 in ADE formation with classical Runge-Kutta scheme
+ ! evolution e1 ! no need since we are just considering shear attenuation
+ if(stage_time_scheme == 4) then
+ e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2 - &
+ e11(i,j,ispec,i_sls) * tauinvnu2)
+ if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
- if(i_stage==1)then
- e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
- endif
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) &
- + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
- else if(i_stage==4)then
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
- (e11_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,2) + &
+ if(i_stage==1) e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
+ else if(i_stage==4)then
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e11_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,2) + &
2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,3) + e11_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
+ endif
-
- ! evolution e13
- if(stage_time_scheme == 1) then
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) + deltat*e13_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
- e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- e13_accel(i,j,ispec,i_sls) = ((dvx_dzl_n(i,j,ispec) + dvz_dxl_n(i,j,ispec)) * phinu2- &
- e13_veloc(i,j,ispec,i_sls)*tauinvnu2) &
- /(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
- e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- endif
-
-
- if(stage_time_scheme == 6) then
- e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) &
- + deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) &
- - deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)
- endif
-
- if(stage_time_scheme == 4) then
- e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2- &
- e13(i,j,ispec,i_sls) * tauinvnu2)
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
- if(i_stage==1)then
- e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
- endif
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) &
- + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
- else if(i_stage==4)then
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
- (e13_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,2) + &
- 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
-
- enddo
- enddo
- enddo
+ e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2 - &
+ e13(i,j,ispec,i_sls) * tauinvnu2)
+ if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
+ if(i_stage==1) e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
+ else if(i_stage==4)then
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e13_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,2) + &
+ 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
+ endif
+ endif
+ enddo
+ enddo; enddo
endif
enddo
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90 2013-09-30 17:22:41 UTC (rev 22905)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90 2013-10-01 01:02:10 UTC (rev 22906)
@@ -46,16 +46,13 @@
source_type,it,NSTEP,anyabs, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
- accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displw_poroelastic,&
- b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
+ accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displs_poroelastic_old,&
+ displw_poroelastic,b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
- jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
- e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu2,&
- phi_nu2,Mu_nu2,N_SLS, &
- rx_viscous,rz_viscous,theta_e,theta_s,&
- b_viscodampx,b_viscodampz,&
+ jacobian,source_time_function,sourcearray,adj_sourcearrays, &
+ e11,e13,hprime_xx,hprimewgll_xx,hprime_zz,hprimewgll_zz,wxgll,wzgll,&
+ inv_tau_sigma_nu2,phi_nu2,Mu_nu2,N_SLS, &
+ rx_viscous,rz_viscous,theta_e,theta_s,b_viscodampx,b_viscodampz,&
ibegin_bottom_poro,iend_bottom_poro,ibegin_top_poro,iend_top_poro, &
ibegin_left_poro,iend_left_poro,ibegin_right_poro,iend_right_poro,&
mufr_k,B_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
@@ -95,7 +92,8 @@
logical, dimension(nspec) :: poroelastic
logical, dimension(4,nelemabs) :: codeabs
- real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: accels_poroelastic,velocs_poroelastic,displs_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: accels_poroelastic,velocs_poroelastic,&
+ displs_poroelastic,displs_poroelastic_old
real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: velocw_poroelastic,displw_poroelastic
real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic
double precision, dimension(2,numat) :: density
@@ -115,8 +113,6 @@
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11,e13
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_veloc,e13_veloc
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_accel,e13_accel
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_LDDRK,e13_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS,stage_time_scheme) :: e11_force_RK, e13_force_RK
@@ -125,8 +121,9 @@
real(kind=CUSTOM_REAL) :: e11_sum,e13_sum
integer :: i_sls
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
- dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) ::dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ !nsub1 denote discrete time step n-1
+ dux_dxl_nsub1,duz_dzl_nsub1,duz_dxl_nsub1,dux_dzl_nsub1
! viscous attenuation (poroelastic media)
double precision, dimension(NGLLX,NGLLZ,nspec) :: rx_viscous
@@ -197,132 +194,107 @@
real(kind=CUSTOM_REAL) :: cpIsquare,cpIIsquare,cssquare,cpIl,cpIIl,csl
! for attenuation
- real(kind=CUSTOM_REAL) :: phinu2,tauinvnu2,theta_n_u,theta_n_v
+ real(kind=CUSTOM_REAL) :: phinu2,tauinvnu2,theta_n_u,theta_nsub1_u
+ real(kind=CUSTOM_REAL) :: bb,coef0,coef1,coef2
! implement attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
! compute Grad(displs_poroelastic) at time step n for attenuation
- if(ATTENUATION_VISCOELASTIC_SOLID) call compute_gradient_attenuation(displs_poroelastic,dux_dxl_n,duz_dxl_n, &
- dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,nglob)
+ call compute_gradient_attenuation(displs_poroelastic,dux_dxl_n,duz_dxl_n, &
+ dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,nglob)
-! compute Grad(velocs_poroelastic) at time step n for attenuation
- call compute_gradient_attenuation(velocs_poroelastic,dvx_dxl_n,dvz_dxl_n, &
- dvx_dzl_n,dvz_dzl_n,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,nglob)
+! compute Grad(displs_poroelastic) at time step n-1 for attenuation
+ call compute_gradient_attenuation(displs_poroelastic_old,dux_dxl_nsub1,duz_dxl_nsub1, &
+ dux_dzl_nsub1,duz_dzl_nsub1,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,nglob)
! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
! loop over spectral elements
do ispec = 1,nspec
if (poroelastic(ispec)) then
-
- do j=1,NGLLZ
- do i=1,NGLLX
+ do j=1,NGLLZ; do i=1,NGLLX
theta_n_u = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
- theta_n_v = dvx_dxl_n(i,j,ispec) + dvz_dzl_n(i,j,ispec)
+ theta_nsub1_u = dux_dxl_nsub1(i,j,ispec) + duz_dzl_nsub1(i,j,ispec)
-! loop on all the standard linear solids
+ ! loop on all the standard linear solids
do i_sls = 1,N_SLS
+ phinu2 = phi_nu2(i,j,ispec,i_sls)
+ tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
+ ! update e1, e11, e13 in convolution formation with modified recursive convolution scheme on basis of
+ ! second-order accurate convolution term calculation from equation (21) of
+ ! Shumin Wang, Robert Lee, and Fernando L. Teixeira,
+ ! Anisotropic-Medium PML for Vector FETD With Modified Basis Functions,
+ ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006)
+ ! evolution e1 ! no need since we are just considering shear attenuation
+ if(stage_time_scheme == 1) then
+ bb = tauinvnu2; coef0 = exp(-bb * deltat)
+ if ( abs(bb) > 1e-5_CUSTOM_REAL ) then
+ coef1 = (1._CUSTOM_REAL - exp(-bb * deltat / 2._CUSTOM_REAL)) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-bb* deltat / 2._CUSTOM_REAL)) * exp(-bb * deltat / 2._CUSTOM_REAL)/ bb
+ else
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
+ endif
-! evolution e1 ! no need since we are just considering shear attenuation
-! if(stage_time_scheme == 1) then
-! Un = e1(i,j,ispec,i_sls)
-! tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
-! tauinvsquare = tauinv * tauinv
-! tauinvcube = tauinvsquare * tauinv
-! tauinvUn = tauinv * Un
-! Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i,j,ispec,i_sls)
-! Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i,j,ispec,i_sls)
-! Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
-! twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
-! fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
-! deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
-! e1(i,j,ispec,i_sls) = Unp1
-! endif
+ e11(i,j,ispec,i_sls) = coef0 * e11(i,j,ispec,i_sls) + &
+ phinu2 * (coef1 * (dux_dxl_n(i,j,ispec)-theta_n_u/TWO) + &
+ coef2 * (dux_dxl_nsub1(i,j,ispec)-theta_nsub1_u/TWO))
-! evolution e11
- if(stage_time_scheme == 1) then
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) + deltat*e11_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e11_accel(i,j,ispec,i_sls)
- e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- e11_accel(i,j,ispec,i_sls) = ((dvx_dxl_n(i,j,ispec)-theta_n_v/TWO) * phinu2- &
- e11_veloc(i,j,ispec,i_sls)*tauinvnu2) &
- /(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
- e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- endif
+ e13(i,j,ispec,i_sls) = coef0 * e13(i,j,ispec,i_sls) + &
+ phinu2 * (coef1 * (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) + &
+ coef2 * (dux_dzl_nsub1(i,j,ispec) + duz_dxl_nsub1(i,j,ispec)))
+ endif
- if(stage_time_scheme == 6) then
- e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
- + deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) &
- - deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*e11_LDDRK(i,j,ispec,i_sls)
- endif
+ ! update e1, e11, e13 in ADE formation with LDDRK scheme
+ ! evolution e1 ! no need since we are just considering shear attenuation
+ if(stage_time_scheme == 6) then
+ e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) + &
+ deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) - &
+ deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
+ e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*e11_LDDRK(i,j,ispec,i_sls)
- if(stage_time_scheme == 4) then
- e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2- &
- e11(i,j,ispec,i_sls) * tauinvnu2)
+ e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) + &
+ deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) - &
+ deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
+ e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)
+ endif
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
- if(i_stage==1)then
- e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
- endif
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) &
- + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
- else if(i_stage==4)then
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
- (e11_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,2) + &
+ ! update e1, e11, e13 in ADE formation with classical Runge-Kutta scheme
+ ! evolution e1 ! no need since we are just considering shear attenuation
+ if(stage_time_scheme == 4) then
+ e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2 - &
+ e11(i,j,ispec,i_sls) * tauinvnu2)
+ if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
+
+ if(i_stage==1) e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
+ else if(i_stage==4)then
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e11_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,2) + &
2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,3) + e11_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
+ endif
-
-! evolution e13
- if(stage_time_scheme == 1) then
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) + deltat*e13_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
- e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- e13_accel(i,j,ispec,i_sls) = ((dvx_dzl_n(i,j,ispec) + dvz_dxl_n(i,j,ispec)) * phinu2- &
- e13_veloc(i,j,ispec,i_sls)*tauinvnu2) &
- /(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
- e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- endif
-
-
- if(stage_time_scheme == 6) then
- e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) &
- + deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) &
- - deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)
- endif
-
- if(stage_time_scheme == 4) then
- e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2- &
- e13(i,j,ispec,i_sls) * tauinvnu2)
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
- if(i_stage==1)then
- e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
- endif
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) &
- + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
- else if(i_stage==4)then
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
- (e13_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,2) + &
- 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
-
- enddo
- enddo
- enddo
+ e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2 - &
+ e13(i,j,ispec,i_sls) * tauinvnu2)
+ if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
+ if(i_stage==1) e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
+ else if(i_stage==4)then
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e13_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,2) + &
+ 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
+ endif
+ endif
+ enddo
+ enddo; enddo
endif
enddo
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-09-30 17:22:41 UTC (rev 22905)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-10-01 01:02:10 UTC (rev 22906)
@@ -47,28 +47,26 @@
source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
- accel_elastic,veloc_elastic,displ_elastic, &
+ accel_elastic,veloc_elastic,displ_elastic,displ_elastic_old, &
density,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,c12ext,c23ext,c25ext,anisotropic,anisotropy, &
- source_time_function,sourcearray,adj_sourcearrays,e1,e11, &
- e13,e1_veloc,e11_veloc,e13_veloc,e1_accel,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
+ source_time_function,sourcearray,adj_sourcearrays, &
+ e1,e11,e13,e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
+ e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, &
+ hprime_xx,hprimewgll_xx,hprime_zz,hprimewgll_zz,wxgll,wzgll, &
+ inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
deltat,coord,add_Bielak_conditions, &
x0_source, z0_source, A_plane, B_plane, C_plane, anglesource_refl, c_inc, c_refl, time_offset,f0, &
v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot,&
nleft,nright,nbot,over_critical_angle,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,b_absorb_elastic_left,&
b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top,nspec_left,nspec_right,&
- nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,&
- e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
- e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, &
+ nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,&
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,nadj_rec_local, &
is_PML,nspec_PML,spec_to_PML,region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
rmemory_displ_elastic,rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
- rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
- rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
+ rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,backward_simulation,STACEY_BOUNDARY_CONDITIONS)
@@ -109,7 +107,7 @@
logical, dimension(nspec) :: elastic,anisotropic
logical, dimension(4,nelemabs) :: codeabs
- real(kind=CUSTOM_REAL), dimension(3,nglob) :: accel_elastic,veloc_elastic,displ_elastic
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: accel_elastic,veloc_elastic,displ_elastic,displ_elastic_old
double precision, dimension(2,numat) :: density
double precision, dimension(4,3,numat) :: poroelastcoef
double precision, dimension(9,numat) :: anisotropy
@@ -128,8 +126,6 @@
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11,e13
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_veloc,e11_veloc,e13_veloc
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_accel,e11_accel,e13_accel
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_LDDRK,e11_LDDRK,e13_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_initial_rk,e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS,stage_time_scheme) :: e1_force_RK, e11_force_RK, e13_force_RK
@@ -137,10 +133,10 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL) :: e1_sum,e11_sum,e13_sum
integer :: i_sls
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) ::dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ !nsub1 denote discrete time step n-1
+ dux_dxl_nsub1,duz_dzl_nsub1,duz_dxl_nsub1,dux_dzl_nsub1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
- dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n
-
! derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
@@ -193,7 +189,7 @@
lambdal_relaxed_viscoelastic,mul_relaxed_viscoelastic,lambdalplusmul_relaxed_viscoel
! for attenuation
- real(kind=CUSTOM_REAL) :: phinu1,phinu2,tauinvnu1,tauinvnu2,theta_n_u,theta_n_v
+ real(kind=CUSTOM_REAL) :: phinu1,phinu2,tauinvnu1,tauinvnu2,theta_n_u,theta_nsub1_u
! for anisotropy
double precision :: c11,c15,c13,c33,c35,c55,c12,c23,c25
@@ -243,162 +239,139 @@
real(kind=CUSTOM_REAL), dimension(3,nglob) :: displ_elastic_new
logical :: backward_simulation
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! implement attenuation
+!!!update momeory variable in viscoelastic simulation
+
if(ATTENUATION_VISCOELASTIC_SOLID) then
- ! compute Grad(displ_elastic) at time step n for attenuation
- call compute_gradient_attenuation(displ_elastic,dux_dxl_n,duz_dxl_n, &
+ ! compute Grad(displ_elastic) at time step n for attenuation
+ call compute_gradient_attenuation(displ_elastic,dux_dxl_n,duz_dxl_n, &
dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
- ! compute Grad(veloc_elastic) at time step n for attenuation
- call compute_gradient_attenuation(veloc_elastic,dvx_dxl_n,dvz_dxl_n, &
- dvx_dzl_n,dvz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
+ ! compute Grad(veloc_elastic) at time step n for attenuation
+ call compute_gradient_attenuation(displ_elastic_old,dux_dxl_nsub1,duz_dxl_nsub1, &
+ dux_dzl_nsub1,duz_dzl_nsub1,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
+ ! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
+ ! loop over spectral elements
+ do ispec = 1,nspec
+ if((.not. PML_BOUNDARY_CONDITIONS) .or. (PML_BOUNDARY_CONDITIONS .and. (.not. is_PML(ispec))))then
+ do j=1,NGLLZ; do i=1,NGLLX
+ theta_n_u = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
+ theta_nsub1_u = dux_dxl_nsub1(i,j,ispec) + duz_dzl_nsub1(i,j,ispec)
- ! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
- ! loop over spectral elements
- do ispec = 1,nspec
+ ! loop on all the standard linear solids
+ do i_sls = 1,N_SLS
+ phinu1 = phi_nu1(i,j,ispec,i_sls)
+ tauinvnu1 = inv_tau_sigma_nu1(i,j,ispec,i_sls)
+ phinu2 = phi_nu2(i,j,ispec,i_sls)
+ tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
- if((.not. PML_BOUNDARY_CONDITIONS) .or. &
- (PML_BOUNDARY_CONDITIONS .and. (.not. is_PML(ispec))))then
+ ! update e1, e11, e13 in convolution formation with modified recursive convolution scheme on basis of
+ ! second-order accurate convolution term calculation from equation (21) of
+ ! Shumin Wang, Robert Lee, and Fernando L. Teixeira,
+ ! Anisotropic-Medium PML for Vector FETD With Modified Basis Functions,
+ ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006)
+ if(stage_time_scheme == 1) then
+ bb = tauinvnu1; coef0 = exp(- bb * deltat)
+ if ( abs(bb) > 1e-5_CUSTOM_REAL ) then
+ coef1 = (1._CUSTOM_REAL - exp(-bb * deltat / 2._CUSTOM_REAL)) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-bb* deltat / 2._CUSTOM_REAL)) * exp(-bb * deltat / 2._CUSTOM_REAL)/ bb
+ else
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
+ endif
- do j=1,NGLLZ
- do i=1,NGLLX
- theta_n_u = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
- theta_n_v = dvx_dxl_n(i,j,ispec) + dvz_dzl_n(i,j,ispec)
+ e1(i,j,ispec,i_sls) = coef0 * e1(i,j,ispec,i_sls) + &
+ phinu1 * (coef1 * theta_n_u + coef2 * theta_nsub1_u)
- ! loop on all the standard linear solids
- do i_sls = 1,N_SLS
+ bb = tauinvnu2; coef0 = exp(-bb * deltat)
+ if ( abs(bb) > 1e-5_CUSTOM_REAL ) then
+ coef1 = (1._CUSTOM_REAL - exp(-bb * deltat / 2._CUSTOM_REAL)) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-bb* deltat / 2._CUSTOM_REAL)) * exp(-bb * deltat / 2._CUSTOM_REAL)/ bb
+ else
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
+ endif
- phinu1 = phi_nu1(i,j,ispec,i_sls)
- tauinvnu1 = inv_tau_sigma_nu1(i,j,ispec,i_sls)
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
+ e11(i,j,ispec,i_sls) = coef0 * e11(i,j,ispec,i_sls) + &
+ phinu2 * (coef1 * (dux_dxl_n(i,j,ispec)-theta_n_u/TWO) + &
+ coef2 * (dux_dxl_nsub1(i,j,ispec)-theta_nsub1_u/TWO))
- ! evolution e1
- if(stage_time_scheme == 1) then
- e1(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls) + deltat*e1_veloc(i,j,ispec,i_sls) + &
- deltatsquareover2*e1_accel(i,j,ispec,i_sls)
- e1_veloc(i,j,ispec,i_sls) = e1_veloc(i,j,ispec,i_sls) + deltatover2*e1_accel(i,j,ispec,i_sls)
- e1_accel(i,j,ispec,i_sls) = (theta_n_v * phinu1 - e1_veloc(i,j,ispec,i_sls) * tauinvnu1) / &
- (1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu1*deltat)
- e1_veloc(i,j,ispec,i_sls) = e1_veloc(i,j,ispec,i_sls) + deltatover2*e1_accel(i,j,ispec,i_sls)
- endif
+ e13(i,j,ispec,i_sls) = coef0 * e13(i,j,ispec,i_sls) + &
+ phinu2 * (coef1 * (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) + &
+ coef2 * (dux_dzl_nsub1(i,j,ispec) + duz_dxl_nsub1(i,j,ispec)))
+ endif
- if(stage_time_scheme == 6) then
- e1_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls) + &
- deltat * (theta_n_u * phinu1 - e1(i,j,ispec,i_sls) * tauinvnu1)
- e1(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls) + beta_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls)
- endif
+ ! update e1, e11, e13 in ADE formation with LDDRK scheme
+ if(stage_time_scheme == 6) then
+ e1_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls) + &
+ deltat * (theta_n_u * phinu1 - e1(i,j,ispec,i_sls) * tauinvnu1)
+ e1(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls) + beta_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls)
- if(stage_time_scheme == 4) then
- e1_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (theta_n_u * phinu1 - e1(i,j,ispec,i_sls) * tauinvnu1)
+ e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) + &
+ deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) - &
+ deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
+ e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*e11_LDDRK(i,j,ispec,i_sls)
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
+ e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) + &
+ deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) - &
+ deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
+ e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)
+ endif
- if(i_stage==1)then
- e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
- endif
+ ! update e1, e11, e13 in ADE formation with classical Runge-Kutta scheme
+ if(stage_time_scheme == 4) then
+ e1_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (theta_n_u * phinu1 - e1(i,j,ispec,i_sls) * tauinvnu1)
- e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) &
- + weight_rk * e1_force_RK(i,j,ispec,i_sls,i_stage)
+ if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
- else if(i_stage==4)then
+ if(i_stage==1) e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
+ e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + weight_rk * e1_force_RK(i,j,ispec,i_sls,i_stage)
+ else if(i_stage==4)then
+ e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e1_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e1_force_RK(i,j,ispec,i_sls,2) + &
+ 2._CUSTOM_REAL * e1_force_RK(i,j,ispec,i_sls,3) + e1_force_RK(i,j,ispec,i_sls,4))
+ endif
- e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
- (e1_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e1_force_RK(i,j,ispec,i_sls,2) + &
- 2._CUSTOM_REAL * e1_force_RK(i,j,ispec,i_sls,3) + e1_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
+ e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2 - &
+ e11(i,j,ispec,i_sls) * tauinvnu2)
+ if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
-
- ! evolution e11
- if(stage_time_scheme == 1) then
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) + deltat*e11_veloc(i,j,ispec,i_sls) + &
- deltatsquareover2*e11_accel(i,j,ispec,i_sls)
- e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- e11_accel(i,j,ispec,i_sls) = ((dvx_dxl_n(i,j,ispec)-theta_n_v/TWO) * phinu2- &
- e11_veloc(i,j,ispec,i_sls)*tauinvnu2) / &
- (1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
- e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- endif
-
- if(stage_time_scheme == 6) then
- e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
- + deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) &
- - deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*e11_LDDRK(i,j,ispec,i_sls)
- endif
-
- if(stage_time_scheme == 4) then
- e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2- &
- e11(i,j,ispec,i_sls) * tauinvnu2)
-
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
- if(i_stage==1)then
- e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
- endif
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) &
- + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
- else if(i_stage==4)then
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
- (e11_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,2) + &
+ if(i_stage==1) e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
+ else if(i_stage==4)then
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e11_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,2) + &
2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,3) + e11_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
+ endif
- ! evolution e13
- if(stage_time_scheme == 1) then
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) + deltat*e13_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
- e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- e13_accel(i,j,ispec,i_sls) = ((dvx_dzl_n(i,j,ispec) + dvz_dxl_n(i,j,ispec)) * phinu2- &
- e13_veloc(i,j,ispec,i_sls)*tauinvnu2) / &
- (1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
- e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- endif
+ e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2 - &
+ e13(i,j,ispec,i_sls) * tauinvnu2)
+ if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
+ if(i_stage==1) e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
+ else if(i_stage==4)then
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e13_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,2) + &
+ 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
+ endif
+ endif
+ enddo
+ enddo; enddo
+ endif
+ enddo
+ endif
+!!!! end of update momeory variable in viscoelastic simulation
-
- if(stage_time_scheme == 6) then
- e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) &
- + deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) &
- - deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)
- endif
-
- if(stage_time_scheme == 4) then
- e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2- &
- e13(i,j,ispec,i_sls) * tauinvnu2)
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
- if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
- if(i_stage==1)then
- e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
- endif
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) &
- + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
- else if(i_stage==4)then
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
- (e13_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,2) + &
- 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
- enddo
- enddo
- enddo
- endif
- enddo
- endif ! end of test on attenuation
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
! this to avoid a warning at execution time about an undefined variable being used
! for the SH component in the case of a P-SV calculation, and vice versa
sigma_xx = 0
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-09-30 17:22:41 UTC (rev 22905)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-10-01 01:02:10 UTC (rev 22906)
@@ -473,7 +473,7 @@
! material properties of the elastic medium
double precision :: mul_unrelaxed_elastic,lambdal_unrelaxed_elastic,lambdaplus2mu_unrelaxed_elastic,kappal
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_elastic,veloc_elastic,displ_elastic
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_elastic,veloc_elastic,displ_elastic,displ_elastic_old
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: veloc_elastic_LDDRK,displ_elastic_LDDRK,&
veloc_elastic_LDDRK_temp
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: accel_elastic_rk,veloc_elastic_rk
@@ -484,7 +484,7 @@
! material properties of the poroelastic medium (solid phase:s and fluid phase [defined as w=phi(u_f-u_s)]: w)
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
- accels_poroelastic,velocs_poroelastic,displs_poroelastic
+ accels_poroelastic,velocs_poroelastic,displs_poroelastic, displs_poroelastic_old
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
velocs_poroelastic_LDDRK,displs_poroelastic_LDDRK
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
@@ -595,10 +595,6 @@
integer nspec_allocate
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
-!DK DK e1_veloc,e11_veloc,e13_veloc denote first derivative of e1,e11,e13 respect to t
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_veloc,e11_veloc,e13_veloc
-!DK DK e1_accel,e11_accel,e13_accel denote second derivative of e1,e11,e13 respect to t
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_accel,e11_accel,e13_accel
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_LDDRK,e11_LDDRK,e13_LDDRK
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_initial_rk,e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: e1_force_rk,e11_force_rk,e13_force_rk
@@ -607,11 +603,6 @@
real(kind=CUSTOM_REAL), dimension(:,:,:) , allocatable :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL) :: Mu_nu1_sent,Mu_nu2_sent
-
-
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
- dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n
-
! for viscous attenuation
double precision, dimension(:,:,:), allocatable :: &
rx_viscous,rz_viscous,viscox,viscoz
@@ -993,6 +984,9 @@
integer :: i_stage,stage_time_scheme
real(kind=CUSTOM_REAL), dimension(Nstages):: alpha_LDDRK,beta_LDDRK,c_LDDRK
+ ! parameters used in LDDRK scheme, from equation (2) of
+ ! Berland, J., Bogey, C., & Bailly, C.
+ ! Low-dissipation and low-dispersion fourth-order Runge–Kutta algorithm, Computers & Fluids, 35(10), 1459-1463.
Data alpha_LDDRK /0.0_CUSTOM_REAL,-0.737101392796_CUSTOM_REAL, &
-1.634740794341_CUSTOM_REAL,-0.744739003780_CUSTOM_REAL, &
-1.469897351522_CUSTOM_REAL,-2.813971388035_CUSTOM_REAL/
@@ -1376,26 +1370,10 @@
allocate(e11(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(e13(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e1_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e11_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e13_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
-
- allocate(e1_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e11_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e13_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
-
e1(:,:,:,:) = 0._CUSTOM_REAL
e11(:,:,:,:) = 0._CUSTOM_REAL
e13(:,:,:,:) = 0._CUSTOM_REAL
- e1_veloc(:,:,:,:) = 0._CUSTOM_REAL
- e11_veloc(:,:,:,:) = 0._CUSTOM_REAL
- e13_veloc(:,:,:,:) = 0._CUSTOM_REAL
-
- e1_accel(:,:,:,:) = 0._CUSTOM_REAL
- e11_accel(:,:,:,:) = 0._CUSTOM_REAL
- e13_accel(:,:,:,:) = 0._CUSTOM_REAL
-
if(time_stepping_scheme == 2)then
allocate(e1_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(e11_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
@@ -1430,15 +1408,6 @@
e1_force_rk(:,:,:,:,:) = 0._CUSTOM_REAL
e11_force_rk(:,:,:,:,:) = 0._CUSTOM_REAL
e13_force_rk(:,:,:,:,:) = 0._CUSTOM_REAL
-
- allocate(dux_dxl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(duz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(duz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dux_dzl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dvx_dxl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dvz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dvz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dvx_dzl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(Mu_nu1(NGLLX,NGLLZ,nspec))
allocate(Mu_nu2(NGLLX,NGLLZ,nspec))
endif
@@ -2659,6 +2628,7 @@
nglob_elastic = 1
endif
allocate(displ_elastic(3,nglob_elastic))
+ allocate(displ_elastic_old(3,nglob_elastic))
allocate(veloc_elastic(3,nglob_elastic))
allocate(accel_elastic(3,nglob_elastic))
allocate(accel_elastic_adj_coupling(3,nglob_elastic))
@@ -2667,16 +2637,16 @@
allocate(rmass_inverse_elastic_three(nglob_elastic))
if(time_stepping_scheme==2)then
- allocate(displ_elastic_LDDRK(3,nglob_elastic))
- allocate(veloc_elastic_LDDRK(3,nglob_elastic))
- allocate(veloc_elastic_LDDRK_temp(3,nglob_elastic))
+ allocate(displ_elastic_LDDRK(3,nglob_elastic))
+ allocate(veloc_elastic_LDDRK(3,nglob_elastic))
+ allocate(veloc_elastic_LDDRK_temp(3,nglob_elastic))
endif
if(time_stepping_scheme == 3)then
- allocate(accel_elastic_rk(3,nglob_elastic,stage_time_scheme))
- allocate(veloc_elastic_rk(3,nglob_elastic,stage_time_scheme))
- allocate(veloc_elastic_initial_rk(3,nglob_elastic))
- allocate(displ_elastic_initial_rk(3,nglob_elastic))
+ allocate(accel_elastic_rk(3,nglob_elastic,stage_time_scheme))
+ allocate(veloc_elastic_rk(3,nglob_elastic,stage_time_scheme))
+ allocate(veloc_elastic_initial_rk(3,nglob_elastic))
+ allocate(displ_elastic_initial_rk(3,nglob_elastic))
endif
! extra array if adjoint and kernels calculation
@@ -2729,6 +2699,7 @@
nglob_poroelastic = 1
endif
allocate(displs_poroelastic(NDIM,nglob_poroelastic))
+ allocate(displs_poroelastic_old(NDIM,nglob_poroelastic))
allocate(velocs_poroelastic(NDIM,nglob_poroelastic))
allocate(accels_poroelastic(NDIM,nglob_poroelastic))
allocate(accels_poroelastic_adj_coupling(NDIM,nglob_poroelastic))
@@ -3602,6 +3573,7 @@
! initialize arrays to zero
displ_elastic = 0._CUSTOM_REAL
+ displ_elastic_old = 0._CUSTOM_REAL
veloc_elastic = 0._CUSTOM_REAL
accel_elastic = 0._CUSTOM_REAL
@@ -3625,6 +3597,7 @@
endif
displs_poroelastic = 0._CUSTOM_REAL
+ displs_poroelastic_old = 0._CUSTOM_REAL
velocs_poroelastic = 0._CUSTOM_REAL
accels_poroelastic = 0._CUSTOM_REAL
displw_poroelastic = 0._CUSTOM_REAL
@@ -3632,24 +3605,24 @@
accelw_poroelastic = 0._CUSTOM_REAL
if(time_stepping_scheme == 2) then
- displs_poroelastic_LDDRK = 0._CUSTOM_REAL
- velocs_poroelastic_LDDRK = 0._CUSTOM_REAL
- displw_poroelastic_LDDRK = 0._CUSTOM_REAL
- velocw_poroelastic_LDDRK = 0._CUSTOM_REAL
+ displs_poroelastic_LDDRK = 0._CUSTOM_REAL
+ velocs_poroelastic_LDDRK = 0._CUSTOM_REAL
+ displw_poroelastic_LDDRK = 0._CUSTOM_REAL
+ velocw_poroelastic_LDDRK = 0._CUSTOM_REAL
endif
if(time_stepping_scheme == 3) then
- accels_poroelastic_rk = 0._CUSTOM_REAL
- velocs_poroelastic_rk = 0._CUSTOM_REAL
+ accels_poroelastic_rk = 0._CUSTOM_REAL
+ velocs_poroelastic_rk = 0._CUSTOM_REAL
- accelw_poroelastic_rk = 0._CUSTOM_REAL
- velocw_poroelastic_rk = 0._CUSTOM_REAL
+ accelw_poroelastic_rk = 0._CUSTOM_REAL
+ velocw_poroelastic_rk = 0._CUSTOM_REAL
- velocs_poroelastic_initial_rk = 0._CUSTOM_REAL
- displs_poroelastic_initial_rk = 0._CUSTOM_REAL
+ velocs_poroelastic_initial_rk = 0._CUSTOM_REAL
+ displs_poroelastic_initial_rk = 0._CUSTOM_REAL
- velocw_poroelastic_initial_rk = 0._CUSTOM_REAL
- displw_poroelastic_initial_rk = 0._CUSTOM_REAL
+ velocw_poroelastic_initial_rk = 0._CUSTOM_REAL
+ displw_poroelastic_initial_rk = 0._CUSTOM_REAL
endif
@@ -4926,12 +4899,14 @@
!! DK DK whose dimension is the product of the two dimensions, the second dimension being equal to 1
do i = 1,3*nglob_elastic !! DK DK here change 3 to NDIM when/if we suppress the 2nd component of the arrays (the SH component)
! do i = 1,NDIM*nglob_elastic !! DK DK this should be the correct size in principle, but not here because of the SH component
+ displ_elastic_old(i,1) = displ_elastic(i,1) + deltatsquareover2 * accel_elastic(i,1)
displ_elastic(i,1) = displ_elastic(i,1) &
+ deltat*veloc_elastic(i,1) &
+ deltatsquareover2*accel_elastic(i,1)
veloc_elastic(i,1) = veloc_elastic(i,1) + deltatover2*accel_elastic(i,1)
enddo
#else
+ displ_elastic_old = displ_elastic + deltatsquareover2 * accel_elastic
displ_elastic = displ_elastic &
+ deltat*veloc_elastic &
+ deltatsquareover2*accel_elastic
@@ -4955,6 +4930,7 @@
if(time_stepping_scheme==1)then
!for the solid
+ displs_poroelastic_old = displs_poroelastic + deltatover2 * accels_poroelastic
displs_poroelastic = displs_poroelastic &
+ deltat*velocs_poroelastic &
+ deltatsquareover2*accels_poroelastic
@@ -5747,14 +5723,14 @@
source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
- accel_elastic,veloc_elastic,displ_elastic, &
+ accel_elastic,veloc_elastic,displ_elastic,displ_elastic_old, &
density,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,c12ext,c23ext,c25ext,anisotropic,anisotropy, &
source_time_function,sourcearray,adj_sourcearrays, &
- e1,e11,e13,e1_veloc,e11_veloc,e13_veloc,e1_accel,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1, &
- phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
+ e1,e11,e13,e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
+ e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
+ hprime_xx,hprimewgll_xx,hprime_zz,hprimewgll_zz,wxgll,wzgll,&
+ inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
deltat,coord,add_Bielak_conditions, x_source(1), z_source(1), &
A_plane, B_plane, C_plane, anglesource_refl, c_inc, c_refl, time_offset, f0(1),&
v0x_left(1,it),v0z_left(1,it),v0x_right(1,it),v0z_right(1,it),v0x_bot(1,it),v0z_bot(1,it), &
@@ -5762,9 +5738,7 @@
count_left,count_right,count_bottom,over_critical_angle, &
NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
b_absorb_elastic_left,b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top, &
- nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
- e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
- e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
+ nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,max(1,nadj_rec_local), &
is_PML,nspec_PML,spec_to_PML,region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
@@ -5804,14 +5778,14 @@
source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
- b_accel_elastic,b_veloc_elastic,b_displ_elastic, &
+ b_accel_elastic,b_veloc_elastic,b_displ_elastic,displ_elastic_old, &
density,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,c12ext,c23ext,c25ext,anisotropic,anisotropy, &
source_time_function,sourcearray,adj_sourcearrays, &
- e1,e11,e13,e1_veloc,e11_veloc,e13_veloc,e1_accel,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1, &
- phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
+ e1,e11,e13,e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
+ e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
+ hprime_xx,hprimewgll_xx,hprime_zz,hprimewgll_zz,wxgll,wzgll, &
+ inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
deltat,coord,add_Bielak_conditions, x_source(1), z_source(1), &
A_plane, B_plane, C_plane, anglesource_refl, c_inc, c_refl, time_offset, f0(1),&
v0x_left(1,it),v0z_left(1,it),v0x_right(1,it),v0z_right(1,it),v0x_bot(1,it),v0z_bot(1,it), &
@@ -5819,9 +5793,7 @@
count_left,count_right,count_bottom,over_critical_angle, &
NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
b_absorb_elastic_left,b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top, &
- nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
- e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
- e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
+ nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,max(1,nadj_rec_local), &
is_PML,nspec_PML,spec_to_PML,region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
@@ -6617,16 +6589,13 @@
source_type,it,NSTEP,anyabs, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
- accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displw_poroelastic,&
- b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
+ accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displs_poroelastic_old,&
+ displw_poroelastic,b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
- jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
- e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu2,&
- phi_nu2,Mu_nu2,N_SLS, &
- rx_viscous,rz_viscous,theta_e,theta_s,&
- b_viscodampx,b_viscodampz,&
+ jacobian,source_time_function,sourcearray,adj_sourcearrays, &
+ e11,e13,hprime_xx,hprimewgll_xx,hprime_zz,hprimewgll_zz,wxgll,wzgll,&
+ inv_tau_sigma_nu2,phi_nu2,Mu_nu2,N_SLS, &
+ rx_viscous,rz_viscous,theta_e,theta_s,b_viscodampx,b_viscodampz,&
ibegin_edge1_poro,iend_edge1_poro,ibegin_edge3_poro,iend_edge3_poro, &
ibegin_edge4_poro,iend_edge4_poro,ibegin_edge2_poro,iend_edge2_poro,&
mufr_k,B_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
@@ -6642,15 +6611,12 @@
initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
accelw_poroelastic,velocw_poroelastic,displw_poroelastic,velocs_poroelastic,displs_poroelastic,&
- b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
+ displs_poroelastic_old,b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
- jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
- e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu2,&
- phi_nu2,Mu_nu2,N_SLS, &
- rx_viscous,rz_viscous,theta_e,theta_s,&
- b_viscodampx,b_viscodampz,&
+ jacobian,source_time_function,sourcearray,adj_sourcearrays, &
+ e11,e13,hprime_xx,hprimewgll_xx,hprime_zz,hprimewgll_zz,wxgll,wzgll,&
+ inv_tau_sigma_nu2,phi_nu2,Mu_nu2,N_SLS, &
+ rx_viscous,rz_viscous,theta_e,theta_s,b_viscodampx,b_viscodampz,&
ibegin_edge1_poro,iend_edge1_poro,ibegin_edge3_poro,iend_edge3_poro, &
ibegin_edge4_poro,iend_edge4_poro,ibegin_edge2_poro,iend_edge2_poro,&
C_k,M_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
More information about the CIG-COMMITS
mailing list