[cig-commits] r19471 - in seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk: EXAMPLES/attenuation src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Wed Jan 25 04:27:02 PST 2012
Author: xie.zhinan
Date: 2012-01-25 04:27:01 -0800 (Wed, 25 Jan 2012)
New Revision: 19471
Modified:
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_viscoelastic.f90
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90
Log:
add classical runge_kutta for visoelastic problem
Modified: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D 2012-01-25 10:13:29 UTC (rev 19470)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D 2012-01-25 12:27:01 UTC (rev 19471)
@@ -23,10 +23,10 @@
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 1000 # total number of time steps
-deltat = 3.0e-3 # duration of a time step
+nt = 500 # total number of time steps
+deltat = 2.2e-3 # duration of a time step
USER_T0 = 0.0d0 # use this t0 as earliest starting time rather than the automatically calculated one
-time_stepping_scheme = 2 # 1 = Newmark (2nd order), 2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta), 3 = classical 4th-order 4-stage Runge-Kutta
+time_stepping_scheme = 3 # 1 = Newmark (2nd order), 2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta), 3 = classical 4th-order 4-stage Runge-Kutta
# source parameters
NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
Modified: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2012-01-25 10:13:29 UTC (rev 19470)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2012-01-25 12:27:01 UTC (rev 19471)
@@ -61,6 +61,7 @@
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,mu_k,kappa_k,&
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK,c_LDDRK, & !xiezhinan
+ e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, & !xiezhinan
stage_time_scheme,i_stage) !xiezhinan
@@ -120,6 +121,8 @@
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_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
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
@@ -140,7 +143,7 @@
double precision, dimension(stage) :: alpha_LDDRK,beta_LDDRK,c_LDDRK
!temp variable
- double precision :: temp_LDDRK,temper_LDDRK
+ double precision :: temp_time_scheme,temper_time_scheme,weight_rk
!---
@@ -999,15 +1002,48 @@
e1(i,j,ispec,i_sls) = Unp1
endif
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
+
tauinv = inv_tau_sigma_nu1(i,j,ispec,i_sls)
Un = e1(i,j,ispec,i_sls)
- temp_LDDRK = phi_nu1(i,j,ispec,i_sls)
+ temp_time_scheme = phi_nu1(i,j,ispec,i_sls)
e1_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls) + &
- deltat * (theta_n * temp_LDDRK - Un * tauinv)
+ deltat * (theta_n * temp_time_scheme - Un * tauinv)
e1(i,j,ispec,i_sls) = Un + beta_LDDRK(i_stage) * 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 = e11(i,j,ispec,i_sls)
@@ -1025,14 +1061,47 @@
endif
if(stage_time_scheme == 6) then
- temp_LDDRK = dux_dxl_n(i,j,ispec)-theta_n/TWO
- temper_LDDRK = phi_nu2(i,j,ispec,i_sls)
+ temp_time_scheme = dux_dxl_n(i,j,ispec)-theta_n/TWO
+ temper_time_scheme = phi_nu2(i,j,ispec,i_sls)
e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
- + deltat * (temp_LDDRK * temper_LDDRK)- &
+ + deltat * (temp_time_scheme * temper_time_scheme)- &
deltat * (e11(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
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
+
+ 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 = e13(i,j,ispec,i_sls)
@@ -1050,14 +1119,47 @@
endif
if(stage_time_scheme == 6) then
- temp_LDDRK=dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)
- temper_LDDRK=phi_nu2(i,j,ispec,i_sls)
+ 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_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)+&
- deltat * (temp_LDDRK*temper_LDDRK)- &
+ deltat * (temp_time_scheme*temper_time_scheme)- &
deltat * (e13(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
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
+
+ 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
Modified: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90 2012-01-25 10:13:29 UTC (rev 19470)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90 2012-01-25 12:27:01 UTC (rev 19471)
@@ -530,6 +530,8 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
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
double precision, dimension(:,:,:,:), allocatable :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
double precision, dimension(:), allocatable :: inv_tau_sigma_nu1_sent,phi_nu1_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent
double precision, dimension(:,:,:) , allocatable :: Mu_nu1,Mu_nu2
@@ -1212,6 +1214,21 @@
e13_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
endif
+ if(time_stepping_scheme == 3)then
+ allocate(e1_initial_rk(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e11_initial_rk(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e13_initial_rk(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e1_force_rk(NGLLX,NGLLZ,nspec_allocate,N_SLS,stage_time_scheme))
+ allocate(e11_force_rk(NGLLX,NGLLZ,nspec_allocate,N_SLS,stage_time_scheme))
+ allocate(e13_force_rk(NGLLX,NGLLZ,nspec_allocate,N_SLS,stage_time_scheme))
+ e1_initial_rk(:,:,:,:) = 0._CUSTOM_REAL
+ e11_initial_rk(:,:,:,:) = 0._CUSTOM_REAL
+ e13_initial_rk(:,:,:,:) = 0._CUSTOM_REAL
+ e1_force_rk(:,:,:,:,:) = 0._CUSTOM_REAL
+ e11_force_rk(:,:,:,:,:) = 0._CUSTOM_REAL
+ e13_force_rk(:,:,:,:,:) = 0._CUSTOM_REAL
+ endif
+
allocate(dux_dxl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(duz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(duz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
@@ -4817,6 +4834,7 @@
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,mu_k,kappa_k, &
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK,c_LDDRK, & !xiezhinan
+ e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, & !xiezhinan
stage_time_scheme,i_stage) !xiezhinan
if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
@@ -5461,10 +5479,7 @@
displ_elastic(1,:) = displ_elastic_initial_rk(1,:) + weight_rk * veloc_elastic_RK(1,:,i_stage)
displ_elastic(2,:) = displ_elastic_initial_rk(2,:) + weight_rk * veloc_elastic_RK(2,:,i_stage)
displ_elastic(3,:) = displ_elastic_initial_rk(3,:) + weight_rk * veloc_elastic_RK(3,:,i_stage)
-
-
-
elseif(i_stage==4)then
accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
@@ -5507,8 +5522,6 @@
endif
-
-
if(SIMULATION_TYPE == 2) then
b_accel_elastic(1,:) = b_accel_elastic(1,:) * rmass_inverse_elastic(:)
b_accel_elastic(2,:) = b_accel_elastic(2,:) * rmass_inverse_elastic(:)
More information about the CIG-COMMITS
mailing list