[cig-commits] r19606 - seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Thu Feb 9 10:20:26 PST 2012
Author: xie.zhinan
Date: 2012-02-09 10:20:25 -0800 (Thu, 09 Feb 2012)
New Revision: 19606
Modified:
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/read_databases.f90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90
Log:
add LDDRK for adjoint viscoelastic inversion
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90 2012-02-09 02:18:29 UTC (rev 19605)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90 2012-02-09 18:20:25 UTC (rev 19606)
@@ -62,7 +62,11 @@
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k,&
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, &
- stage_time_scheme,i_stage)
+ stage_time_scheme,i_stage,&
+ b_veloc_elastic,b_e1,b_e11,b_e13,&
+ b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare,&
+ b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,& !xiezhinan
+ b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1,b_e1_LDDRK,b_e11_LDDRK,b_e13_LDDRK)
! compute forces for the elastic elements
@@ -90,6 +94,7 @@
logical :: SAVE_FORWARD
double precision :: deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
+ double precision :: b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare !xiezhinan
double precision, dimension(NSOURCES) :: angleforce
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
@@ -111,7 +116,10 @@
real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP,stage_time_scheme) :: source_time_function
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
- real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_displ_elastic
+! real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_displ_elastic !old
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_veloc_elastic,b_displ_elastic !xiezhinan
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_temp_displ_elastic !xiezhinan
+
! real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays !old
real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,stage_time_scheme,3,NGLLX,NGLLZ) :: adj_sourcearrays !xiezhinan
real(kind=CUSTOM_REAL), dimension(nglob) :: mu_k,kappa_k
@@ -127,16 +135,22 @@
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) :: b_e1,b_e11,b_e13 !xiezhinan
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) :: b_e1_LDDRK,b_e11_LDDRK,b_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
double precision, dimension(NGLLX,NGLLZ,nspec,N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL) :: e1_sum,e11_sum,e13_sum
+ real(kind=CUSTOM_REAL) :: b_e1_sum,b_e11_sum,b_e13_sum !xiezhinan
integer :: i_sls
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: & !xiezhinan
+ b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n, & !xiezhinan
+ b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1 !xiezhinan
! derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
@@ -182,6 +196,7 @@
! for attenuation
real(kind=CUSTOM_REAL) :: Un,Unp1,tauinv,Sn,Snp1,theta_n,theta_np1,tauinvsquare,tauinvcube,tauinvUn
+ real(kind=CUSTOM_REAL) :: b_theta_n,b_theta_np1
! for anisotropy
double precision :: c11,c15,c13,c33,c35,c55
@@ -211,9 +226,15 @@
! compute Grad(displ_elastic) at time step n for attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
- temp_displ_elastic = displ_elastic + deltat * veloc_elastic + 0.5 * deltat**2 * accel_elastic
- 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)
+ if(SIMULATION_TYPE == 1)then
+ temp_displ_elastic = displ_elastic + deltat * veloc_elastic + 0.5 * deltat**2 * accel_elastic
+ 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)
+ else
+ b_temp_displ_elastic = b_displ_elastic + b_deltat * b_veloc_elastic + 0.5 * b_deltat**2 * b_accel_elastic
+ call compute_gradient_attenuation(b_displ_elastic,b_dux_dxl_n,b_duz_dxl_n, &
+ b_dux_dzl_n,b_duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
+ endif
endif
ifirstelem = 1
@@ -352,6 +373,7 @@
lambdalplus2mul_relaxed_viscoel = lambdal_relaxed_viscoelastic + TWO*mul_relaxed_viscoelastic
! compute the stress using the unrelaxed Lame parameters (Carcione 2007 page 125)
+ if(SIMULATION_TYPE == 1) then
sigma_xx = lambdaplus2mu_unrelaxed_elastic*dux_dxl + lambdal_unrelaxed_elastic*duz_dzl
sigma_xz = mul_unrelaxed_elastic*(duz_dxl + dux_dzl)
sigma_zz = lambdaplus2mu_unrelaxed_elastic*duz_dzl + lambdal_unrelaxed_elastic*dux_dxl
@@ -373,7 +395,33 @@
sigma_xz = sigma_xz + mul_relaxed_viscoelastic * e13_sum
sigma_zz = sigma_zz + (lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic) * e1_sum &
- TWO * mul_relaxed_viscoelastic * e11_sum
+ endif
+ ! compute the stress using the unrelaxed Lame parameters (Carcione 2007 page 125) adjoint
+ if(SIMULATION_TYPE == 2) then
+ b_sigma_xx = lambdaplus2mu_unrelaxed_elastic*b_dux_dxl + lambdal_unrelaxed_elastic*b_duz_dzl
+ b_sigma_xz = mul_unrelaxed_elastic*(b_duz_dxl + b_dux_dzl)
+ b_sigma_zz = lambdaplus2mu_unrelaxed_elastic*b_duz_dzl + lambdal_unrelaxed_elastic*b_dux_dxl
+
+ ! add the memory variables using the relaxed parameters (Carcione 2007 page 125)
+ ! beware: there is a bug in Carcione's equation (2c) of his 1993 paper for sigma_zz, we fixed it in the code below
+ b_e1_sum = 0._CUSTOM_REAL
+ b_e11_sum = 0._CUSTOM_REAL
+ b_e13_sum = 0._CUSTOM_REAL
+
+ do i_sls = 1,N_SLS
+ b_e1_sum = b_e1_sum + b_e1(i,j,ispec,i_sls)
+ b_e11_sum = b_e11_sum + b_e11(i,j,ispec,i_sls)
+ b_e13_sum = b_e13_sum + b_e13(i,j,ispec,i_sls)
+ enddo
+
+ b_sigma_xx = b_sigma_xx + (lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic) * b_e1_sum &
+ + TWO * mul_relaxed_viscoelastic * b_e11_sum
+ b_sigma_xz = b_sigma_xz + mul_relaxed_viscoelastic * b_e13_sum
+ b_sigma_zz = b_sigma_zz + (lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic) * b_e1_sum &
+ - TWO * mul_relaxed_viscoelastic * b_e11_sum
+ endif
+
else
! no attenuation
@@ -1040,6 +1088,7 @@
! implement attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
+ if(SIMULATION_TYPE == 1) then
! compute Grad(displ_elastic) at time step n+1 for attenuation
call compute_gradient_attenuation(temp_displ_elastic,dux_dxl_np1,duz_dxl_np1, &
dux_dzl_np1,duz_dzl_np1,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
@@ -1207,7 +1256,178 @@
enddo
enddo
enddo
+ endif
+ if(SIMULATION_TYPE == 2) then
+ ! compute Grad(displ_elastic) at time step n+1 for attenuation
+ call compute_gradient_attenuation(b_temp_displ_elastic,b_dux_dxl_np1,b_duz_dxl_np1, &
+ b_dux_dzl_np1,b_duz_dzl_np1,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
+
+ do j=1,NGLLZ
+ do i=1,NGLLX
+
+ b_theta_n = b_dux_dxl_n(i,j,ispec) + b_duz_dzl_n(i,j,ispec)
+ b_theta_np1 = b_dux_dxl_np1(i,j,ispec) + b_duz_dzl_np1(i,j,ispec)
+
+ ! loop on all the standard linear solids
+ do i_sls = 1,N_SLS
+
+ ! evolution e1
+ if(stage_time_scheme == 1) then
+ Un = b_e1(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu1(i,j,ispec,i_sls)
+ tauinvsquare = tauinv * tauinv
+ tauinvcube = tauinvsquare * tauinv
+ tauinvUn = tauinv * Un
+ Sn = b_theta_n * phi_nu1(i,j,ispec,i_sls)
+ Snp1 = b_theta_np1 * phi_nu1(i,j,ispec,i_sls)
+ Unp1 = Un + (b_deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ b_twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ b_fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ b_deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ b_e1(i,j,ispec,i_sls) = Unp1
+ endif
+
+ if(stage_time_scheme == 6) then
+
+ tauinv = inv_tau_sigma_nu1(i,j,ispec,i_sls)
+ Un = b_e1(i,j,ispec,i_sls)
+ temp_time_scheme = phi_nu1(i,j,ispec,i_sls)
+ b_e1_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * b_e1_LDDRK(i,j,ispec,i_sls) + &
+ b_deltat * (b_theta_n * temp_time_scheme - Un * tauinv)
+ b_e1(i,j,ispec,i_sls) = Un + beta_LDDRK(i_stage) * b_e1_LDDRK(i,j,ispec,i_sls)
+ endif
+
+ if(stage_time_scheme == 4) then
+
+ tauinv = inv_tau_sigma_nu1(i,j,ispec,i_sls)
+ Un = e1(i,j,ispec,i_sls)
+ temp_time_scheme = phi_nu1(i,j,ispec,i_sls)
+ e1_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (theta_n * temp_time_scheme - Un * tauinv)
+
+ if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ if(i_stage == 1)weight_rk = 0.5d0
+ if(i_stage == 2)weight_rk = 0.5d0
+ if(i_stage == 3)weight_rk = 1.0d0
+
+ if(i_stage==1)then
+ e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
+ endif
+
+ 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)
+
+ elseif(i_stage==4)then
+
+ e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
+ (e1_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e1_force_RK(i,j,ispec,i_sls,2) + &
+ 2.0d0 * e1_force_RK(i,j,ispec,i_sls,3) + e1_force_RK(i,j,ispec,i_sls,4))
+ endif
+ endif
+
+
+ ! evolution e11
+ if(stage_time_scheme == 1) then
+ Un = b_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 = (b_dux_dxl_n(i,j,ispec) - b_theta_n/TWO) * phi_nu2(i,j,ispec,i_sls)
+ Snp1 = (b_dux_dxl_np1(i,j,ispec) - b_theta_np1/TWO) * phi_nu2(i,j,ispec,i_sls)
+ Unp1 = Un + (b_deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ b_twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ b_fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ b_deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ b_e11(i,j,ispec,i_sls) = Unp1
+ endif
+
+ if(stage_time_scheme == 6) then
+ temp_time_scheme = b_dux_dxl_n(i,j,ispec)- b_theta_n/TWO
+ temper_time_scheme = phi_nu2(i,j,ispec,i_sls)
+ b_e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * b_e11_LDDRK(i,j,ispec,i_sls) &
+ + b_deltat * (temp_time_scheme * temper_time_scheme)- &
+ b_deltat * (b_e11(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
+ b_e11(i,j,ispec,i_sls) = b_e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*b_e11_LDDRK(i,j,ispec,i_sls)
+ endif
+
+ if(stage_time_scheme == 4) then
+
+ temp_time_scheme = dux_dxl_n(i,j,ispec)-theta_n/TWO
+ temper_time_scheme = phi_nu2(i,j,ispec,i_sls)
+ e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (temp_time_scheme * temper_time_scheme- &
+ e11(i,j,ispec,i_sls) * inv_tau_sigma_nu2(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.5d0
+ if(i_stage == 2)weight_rk = 0.5d0
+ if(i_stage == 3)weight_rk = 1.0d0
+ 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)
+ elseif(i_stage==4)then
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
+ (e11_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e11_force_RK(i,j,ispec,i_sls,2) + &
+ 2.0d0 * e11_force_RK(i,j,ispec,i_sls,3) + e11_force_RK(i,j,ispec,i_sls,4))
+ endif
+ endif
+
+ ! evolution e13
+ if(stage_time_scheme == 1) then
+ Un = b_e13(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 = (b_dux_dzl_n(i,j,ispec) + b_duz_dxl_n(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
+ Snp1 = (b_dux_dzl_np1(i,j,ispec) + b_duz_dxl_np1(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
+ Unp1 = Un + (b_deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ b_twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ b_fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ b_deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ b_e13(i,j,ispec,i_sls) = Unp1
+ endif
+
+ if(stage_time_scheme == 6) then
+ temp_time_scheme=b_dux_dzl_n(i,j,ispec) + b_duz_dxl_n(i,j,ispec)
+ temper_time_scheme=phi_nu2(i,j,ispec,i_sls)
+ b_e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * b_e13_LDDRK(i,j,ispec,i_sls)+&
+ b_deltat * (temp_time_scheme*temper_time_scheme)- &
+ b_deltat * (b_e13(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
+ b_e13(i,j,ispec,i_sls) = b_e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * b_e13_LDDRK(i,j,ispec,i_sls)
+ endif
+
+ if(stage_time_scheme == 4) then
+ temp_time_scheme=dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)
+ temper_time_scheme=phi_nu2(i,j,ispec,i_sls)
+ e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (temp_time_scheme*temper_time_scheme- &
+ e13(i,j,ispec,i_sls) * inv_tau_sigma_nu2(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.5d0
+ if(i_stage == 2)weight_rk = 0.5d0
+ if(i_stage == 3)weight_rk = 1.0d0
+ 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)
+ elseif(i_stage==4)then
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
+ (e13_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e13_force_RK(i,j,ispec,i_sls,2) + &
+ 2.0d0 * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
+ endif
+ endif
+
+ enddo
+
+ enddo
+ enddo
+ enddo
+ endif
+
endif ! end of test on attenuation
end subroutine compute_forces_viscoelastic
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/read_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/read_databases.f90 2012-02-09 02:18:29 UTC (rev 19605)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/read_databases.f90 2012-02-09 18:20:25 UTC (rev 19606)
@@ -264,9 +264,9 @@
if( SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. &
(ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART) ) then
- print*, '*************** WARNING ***************'
- print*, 'Anisotropy & Attenuation & Viscous damping are not presently implemented for adjoint calculations'
- stop
+ ! print*, '*************** WARNING ***************' !xiezhinan
+ ! print*, 'Anisotropy & Attenuation & Viscous damping are not presently implemented for adjoint calculations'
+ ! stop
endif
NTSTEP_BETWEEN_OUTPUT_SEISMO = min(NSTEP,NTSTEP_BETWEEN_OUTPUT_INFO)
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90 2012-02-09 02:18:29 UTC (rev 19605)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90 2012-02-09 18:20:25 UTC (rev 19606)
@@ -543,9 +543,11 @@
double precision :: f0_attenuation
integer nspec_allocate
double precision :: deltatsquare,deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
-
+ double precision :: b_deltatsquare,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_e1,b_e11,b_e13
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_LDDRK,e11_LDDRK,e13_LDDRK
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_e1_LDDRK,b_e11_LDDRK,b_e13_LDDRK !xiezhinan
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
double precision, dimension(:,:,:,:), allocatable :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
@@ -556,6 +558,10 @@
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
+ b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,&
+ b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1
+
! for viscous attenuation
double precision, dimension(:,:,:), allocatable :: &
rx_viscous,rz_viscous,viscox,viscoz
@@ -1251,19 +1257,42 @@
e1(:,:,:,:) = 0._CUSTOM_REAL
e11(:,:,:,:) = 0._CUSTOM_REAL
e13(:,:,:,:) = 0._CUSTOM_REAL
+ if(SIMULATION_TYPE == 2)then
+ allocate(b_e1(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(b_e11(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(b_e13(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ b_e1(:,:,:,:) = 0._CUSTOM_REAL
+ b_e11(:,:,:,:) = 0._CUSTOM_REAL
+ b_e13(:,:,:,:) = 0._CUSTOM_REAL
+ endif
if(time_stepping_scheme == 2)then
allocate(e1_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(e11_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(e13_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ if(SIMULATION_TYPE == 2)then
+ allocate(b_e1_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(b_e11_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(b_e13_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ endif
else
allocate(e1_LDDRK(1,1,1,1))
allocate(e11_LDDRK(1,1,1,1))
allocate(e13_LDDRK(1,1,1,1))
+ if(SIMULATION_TYPE == 2)then
+ allocate(b_e1_LDDRK(1,1,1,1))
+ allocate(b_e11_LDDRK(1,1,1,1))
+ allocate(b_e13_LDDRK(1,1,1,1))
endif
+ endif
e1_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
e11_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
e13_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
+ if(SIMULATION_TYPE == 2)then
+ b_e1_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
+ b_e11_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
+ b_e13_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
+ endif
if(time_stepping_scheme == 3)then
allocate(e1_initial_rk(NGLLX,NGLLZ,nspec_allocate,N_SLS))
@@ -1295,6 +1324,17 @@
allocate(duz_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(duz_dxl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(dux_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
+ if(SIMULATION_TYPE == 2)then
+ allocate(b_dux_dxl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(b_duz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(b_duz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(b_dux_dzl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(b_dux_dxl_np1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(b_duz_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(b_duz_dxl_np1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(b_dux_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
+ endif
+
allocate(Mu_nu1(NGLLX,NGLLZ,nspec))
allocate(Mu_nu2(NGLLX,NGLLZ,nspec))
endif
@@ -3434,13 +3474,24 @@
endif ! initialfield
+ if(SIMULATION_TYPE == 1)then
deltatsquare = deltat * deltat
deltatcube = deltatsquare * deltat
deltatfourth = deltatsquare * deltatsquare
twelvedeltat = 12.d0 * deltat
fourdeltatsquare = 4.d0 * deltatsquare
+ endif
+ if(SIMULATION_TYPE == 2)then
+ b_deltatsquare = b_deltat * b_deltat
+ b_deltatcube = b_deltatsquare * b_deltat
+ b_deltatfourth = b_deltatsquare * b_deltatsquare
+
+ b_twelvedeltat = 12.d0 * b_deltat
+ b_fourdeltatsquare = 4.d0 * b_deltatsquare
+ endif
+
! compute the source time function and store it in a text file
if(.not. initialfield) then
@@ -5187,7 +5238,11 @@
nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k, &
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, &
- stage_time_scheme,i_stage)
+ stage_time_scheme,i_stage,&
+ b_veloc_elastic,b_e1,b_e11,b_e13,&
+ b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare,&
+ b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,& !xiezhinan
+ b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1,b_e1_LDDRK,b_e11_LDDRK,b_e13_LDDRK)
if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
!--- left absorbing boundary
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90 2012-02-09 02:18:29 UTC (rev 19605)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90 2012-02-09 18:20:25 UTC (rev 19606)
@@ -477,12 +477,15 @@
else
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
- call MPI_SEND(sisux(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+! call MPI_SEND(sisux(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror) !old
+ call MPI_SEND(sisux(1,1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror) !xiezhinan
if ( number_of_components >= 2 ) then
- call MPI_SEND(sisuz(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+! call MPI_SEND(sisuz(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror) !old
+ call MPI_SEND(sisuz(1,1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)!xiezhinan
end if
if ( number_of_components == 3 ) then
- call MPI_SEND(siscurl(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+! call MPI_SEND(siscurl(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror) !old
+ call MPI_SEND(siscurl(1,1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)!xiezhinan
end if
end if
#endif
More information about the CIG-COMMITS
mailing list