[cig-commits] r19427 - 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
Tue Jan 24 07:48:11 PST 2012
Author: xie.zhinan
Date: 2012-01-24 07:48:10 -0800 (Tue, 24 Jan 2012)
New Revision: 19427
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/specfem2D.F90
Log:
add runge kutta implementation for elastic case
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-24 12:02:12 UTC (rev 19426)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D 2012-01-24 15:48:10 UTC (rev 19427)
@@ -23,10 +23,10 @@
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 1500 # total number of time steps
-deltat = 7.5e-4 # duration of a time step
+nt = 2000 # total number of time steps
+deltat = 0.75e-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 = 1 # 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 = 2 # 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/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90 2012-01-24 12:02:12 UTC (rev 19426)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90 2012-01-24 15:48:10 UTC (rev 19427)
@@ -436,6 +436,8 @@
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_elastic,veloc_elastic,displ_elastic
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: veloc_elastic_LDDRK,displ_elastic_LDDRK !xiezhinan
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: accel_elastic_RK,veloc_elastic_RK !xiezhinan
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: veloc_elastic_initial_rk,displ_elastic_initial_rk !xiezhinan
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_elastic_adj_coupling
double precision, dimension(:,:), allocatable :: &
coord, flagrange,xinterp,zinterp,Uxinterp,Uzinterp,vector_field_display
@@ -893,6 +895,9 @@
integer :: i_stage,stage_time_scheme
double precision, dimension(6):: alpha_LDDRK,beta_LDDRK,c_LDDRK
+! for rk44
+ double precision :: weight_rk
+
!***********************************************************************
!
! i n i t i a l i z a t i o n p h a s e
@@ -2362,6 +2367,13 @@
allocate(veloc_elastic_LDDRK(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))
+ endif
+
! extra array if adjoint and kernels calculation
if(SIMULATION_TYPE == 2 .and. any_elastic) then
allocate(b_displ_elastic(3,nglob))
@@ -2942,6 +2954,13 @@
veloc_elastic_LDDRK = 0._CUSTOM_REAL
endif
+ if(time_stepping_scheme == 3)then
+ accel_elastic_RK = 0._CUSTOM_REAL
+ veloc_elastic_RK = 0._CUSTOM_REAL
+ veloc_elastic_initial_rk = 0._CUSTOM_REAL
+ displ_elastic_initial_rk = 0._CUSTOM_REAL
+ endif
+
displs_poroelastic = 0._CUSTOM_REAL
velocs_poroelastic = 0._CUSTOM_REAL
accels_poroelastic = 0._CUSTOM_REAL
@@ -4162,6 +4181,7 @@
! update displacement using finite-difference time scheme (Newmark)
if(any_elastic) then
+
if(time_stepping_scheme==1)then
displ_elastic = displ_elastic &
+ deltat*veloc_elastic &
@@ -4170,9 +4190,11 @@
accel_elastic_adj_coupling = accel_elastic
accel_elastic = ZERO
endif
+
if(time_stepping_scheme==2)then
accel_elastic = ZERO
endif
+
if(time_stepping_scheme==3)then
accel_elastic = ZERO
endif
@@ -5373,7 +5395,7 @@
accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic
- veloc_elastic = veloc_elastic + deltatover2*accel_elastic
+ veloc_elastic = veloc_elastic + deltatover2 * accel_elastic
endif
if(time_stepping_scheme == 2)then
@@ -5381,15 +5403,97 @@
accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic
- veloc_elastic_LDDRK=alpha_LDDRK(i_stage)*veloc_elastic_LDDRK+deltat*accel_elastic
- displ_elastic_LDDRK=alpha_LDDRK(i_stage)*displ_elastic_LDDRK+deltat*veloc_elastic
+ veloc_elastic_LDDRK = alpha_LDDRK(i_stage) * veloc_elastic_LDDRK + deltat * accel_elastic
+ displ_elastic_LDDRK = alpha_LDDRK(i_stage) * displ_elastic_LDDRK + deltat * veloc_elastic
- veloc_elastic=veloc_elastic+beta_LDDRK(i_stage)*veloc_elastic_LDDRK
- displ_elastic=displ_elastic+beta_LDDRK(i_stage)*displ_elastic_LDDRK
+ veloc_elastic = veloc_elastic + beta_LDDRK(i_stage) * veloc_elastic_LDDRK
+ displ_elastic = displ_elastic + beta_LDDRK(i_stage) * displ_elastic_LDDRK
+
endif
if(time_stepping_scheme == 3)then
+
+ 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
+ accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
+ accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
+ accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic
+
+ accel_elastic_RK(1,:,i_stage) = deltat * accel_elastic(1,:)
+ accel_elastic_RK(2,:,i_stage) = deltat * accel_elastic(2,:)
+ accel_elastic_RK(3,:,i_stage) = deltat * accel_elastic(3,:)
+
+ veloc_elastic_RK(1,:,i_stage) = deltat * veloc_elastic(1,:)
+ veloc_elastic_RK(2,:,i_stage) = deltat * veloc_elastic(2,:)
+ veloc_elastic_RK(3,:,i_stage) = deltat * veloc_elastic(3,:)
+
+ if(i_stage==1)then
+
+ veloc_elastic_initial_rk(1,:) = veloc_elastic(1,:)
+ veloc_elastic_initial_rk(2,:) = veloc_elastic(2,:)
+ veloc_elastic_initial_rk(3,:) = veloc_elastic(3,:)
+
+ displ_elastic_initial_rk(1,:) = displ_elastic(1,:)
+ displ_elastic_initial_rk(2,:) = displ_elastic(2,:)
+ displ_elastic_initial_rk(3,:) = displ_elastic(3,:)
+
+ endif
+
+ veloc_elastic(1,:) = veloc_elastic_initial_rk(1,:) + weight_rk * accel_elastic_RK(1,:,i_stage)
+ veloc_elastic(2,:) = veloc_elastic_initial_rk(2,:) + weight_rk * accel_elastic_RK(2,:,i_stage)
+ veloc_elastic(3,:) = veloc_elastic_initial_rk(3,:) + weight_rk * accel_elastic_RK(3,:,i_stage)
+
+ 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
+ accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
+ accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic
+
+ accel_elastic_RK(1,:,i_stage) = deltat * accel_elastic(1,:)
+ accel_elastic_RK(2,:,i_stage) = deltat * accel_elastic(2,:)
+ accel_elastic_RK(3,:,i_stage) = deltat * accel_elastic(3,:)
+
+ veloc_elastic_RK(1,:,i_stage) = deltat * veloc_elastic(1,:)
+ veloc_elastic_RK(2,:,i_stage) = deltat * veloc_elastic(2,:)
+ veloc_elastic_RK(3,:,i_stage) = deltat * veloc_elastic(3,:)
+
+ veloc_elastic(1,:) = veloc_elastic_initial_rk(1,:) + 1.0d0 / 6.0d0 * &
+ (accel_elastic_RK(1,:,1) + 2.0d0 * accel_elastic_RK(1,:,2) + &
+ 2.0d0 * accel_elastic_RK(1,:,3) + accel_elastic_RK(1,:,4))
+
+ veloc_elastic(2,:) = veloc_elastic_initial_rk(2,:) + 1.0d0 / 6.0d0 * &
+ (accel_elastic_RK(2,:,1) + 2.0d0 * accel_elastic_RK(2,:,2) + &
+ 2.0d0 * accel_elastic_RK(2,:,3) + accel_elastic_RK(2,:,4))
+
+ veloc_elastic(3,:) = veloc_elastic_initial_rk(3,:) + 1.0d0 / 6.0d0 * &
+ (accel_elastic_RK(3,:,1) + 2.0d0 * accel_elastic_RK(3,:,2) + &
+ 2.0d0 * accel_elastic_RK(3,:,3) + accel_elastic_RK(3,:,4))
+
+ displ_elastic(1,:) = displ_elastic_initial_rk(1,:) + 1.0d0 / 6.0d0 * &
+ (veloc_elastic_RK(1,:,1) + 2.0d0 * veloc_elastic_RK(1,:,2) + &
+ 2.0d0 * veloc_elastic_RK(1,:,3) + veloc_elastic_RK(1,:,4))
+
+ displ_elastic(2,:) = displ_elastic_initial_rk(2,:) + 1.0d0 / 6.0d0 * &
+ (veloc_elastic_RK(2,:,1) + 2.0d0 * veloc_elastic_RK(2,:,2) + &
+ 2.0d0 * veloc_elastic_RK(2,:,3) + veloc_elastic_RK(2,:,4))
+
+ displ_elastic(3,:) = displ_elastic_initial_rk(3,:) + 1.0d0 / 6.0d0 * &
+ (veloc_elastic_RK(3,:,1) + 2.0d0 * veloc_elastic_RK(3,:,2) + &
+ 2.0d0 * veloc_elastic_RK(3,:,3) + veloc_elastic_RK(3,:,4))
+
+ endif
+
endif
More information about the CIG-COMMITS
mailing list