[cig-commits] r19473 - seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Wed Jan 25 07:42:43 PST 2012
Author: xie.zhinan
Date: 2012-01-25 07:42:42 -0800 (Wed, 25 Jan 2012)
New Revision: 19473
Modified:
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90
Log:
add runge kutta and LDDRK time implementation for acoustic part
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 13:41:40 UTC (rev 19472)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90 2012-01-25 15:42:42 UTC (rev 19473)
@@ -436,7 +436,7 @@
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 :: 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 :: &
@@ -464,6 +464,9 @@
! for acoustic medium
real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_dot_acoustic_LDDRK, potential_acoustic_LDDRK
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_acoustic_initial_rk, potential_dot_acoustic_initial_rk
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: potential_dot_dot_acoustic_rk, potential_dot_acoustic_rk
real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_acoustic_adj_coupling
! inverse mass matrices
@@ -2395,8 +2398,8 @@
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(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
@@ -2576,7 +2579,19 @@
allocate(potential_dot_acoustic(nglob_acoustic))
allocate(potential_dot_dot_acoustic(nglob_acoustic))
allocate(rmass_inverse_acoustic(nglob_acoustic))
+ if(time_stepping_scheme == 2) then
+ allocate(potential_acoustic_LDDRK(nglob_acoustic))
+ allocate(potential_dot_acoustic_LDDRK(nglob_acoustic))
+ endif
+ if(time_stepping_scheme == 3) then
+ allocate(potential_acoustic_initial_rk(nglob_acoustic))
+ allocate(potential_dot_acoustic_initial_rk(nglob_acoustic))
+ allocate(potential_dot_dot_acoustic_rk(nglob_acoustic,stage_time_scheme))
+ allocate(potential_dot_acoustic_rk(nglob_acoustic,stage_time_scheme))
+ endif
+
+
if(SIMULATION_TYPE == 2 .and. any_acoustic) then
allocate(b_potential_acoustic(nglob))
allocate(b_potential_dot_acoustic(nglob))
@@ -2982,8 +2997,8 @@
endif
if(time_stepping_scheme == 3)then
- accel_elastic_RK = 0._CUSTOM_REAL
- veloc_elastic_RK = 0._CUSTOM_REAL
+ 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
@@ -2999,6 +3014,18 @@
potential_dot_acoustic = 0._CUSTOM_REAL
potential_dot_dot_acoustic = 0._CUSTOM_REAL
+ if(time_stepping_scheme == 2 )then
+ potential_acoustic_LDDRK = 0._CUSTOM_REAL
+ potential_dot_acoustic_LDDRK = 0._CUSTOM_REAL
+ endif
+
+ if(time_stepping_scheme == 3 )then
+ potential_acoustic_initial_rk = 0._CUSTOM_REAL
+ potential_dot_acoustic_initial_rk = 0._CUSTOM_REAL
+ potential_dot_dot_acoustic_rk = 0._CUSTOM_REAL
+ potential_dot_acoustic_rk = 0._CUSTOM_REAL
+ endif
+
!
!----- Files where viscous damping are saved during forward wavefield calculation
!
@@ -4775,10 +4802,64 @@
! ************************************************************************************
if(any_acoustic) then
-
+ if(time_stepping_scheme == 1)then
potential_dot_dot_acoustic = potential_dot_dot_acoustic * rmass_inverse_acoustic
potential_dot_acoustic = potential_dot_acoustic + deltatover2*potential_dot_dot_acoustic
+ endif
+ if(time_stepping_scheme == 2)then
+
+ potential_dot_dot_acoustic = potential_dot_dot_acoustic * rmass_inverse_acoustic
+
+ potential_dot_acoustic_LDDRK = alpha_LDDRK(i_stage) * potential_dot_acoustic_LDDRK &
+ + deltat * potential_dot_dot_acoustic
+
+ potential_acoustic_LDDRK = alpha_LDDRK(i_stage) * potential_acoustic_LDDRK &
+ +deltat*potential_dot_acoustic
+
+ potential_dot_acoustic = potential_dot_acoustic + beta_LDDRK(i_stage) * potential_dot_acoustic_LDDRK
+
+ potential_acoustic = potential_acoustic + beta_LDDRK(i_stage) * potential_acoustic_LDDRK
+
+ endif
+
+ if(time_stepping_scheme == 3)then
+
+ potential_dot_dot_acoustic = potential_dot_dot_acoustic * rmass_inverse_acoustic
+
+ potential_dot_dot_acoustic_rk(:,i_stage) = deltat * potential_dot_dot_acoustic(:)
+ potential_dot_acoustic_rk(:,i_stage) = deltat * potential_dot_acoustic(:)
+
+ 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
+
+ potential_dot_acoustic_initial_rk(:) = potential_dot_acoustic(:)
+ potential_acoustic_initial_rk(:) = potential_acoustic(:)
+
+ endif
+
+ potential_dot_acoustic(:) = potential_dot_acoustic_initial_rk(:) + weight_rk * potential_dot_dot_acoustic_rk(:,i_stage)
+ potential_acoustic(:) = potential_acoustic_initial_rk(:) + weight_rk * potential_dot_acoustic_rk(:,i_stage)
+
+ elseif(i_stage==4)then
+
+ potential_dot_acoustic(:) = potential_dot_acoustic_initial_rk(:) + 1.0d0 / 6.0d0 * &
+ (potential_dot_dot_acoustic_rk(:,1) + 2.0d0 * potential_dot_dot_acoustic_rk(:,2) + &
+ 2.0d0 * potential_dot_dot_acoustic_rk(:,3) + potential_dot_dot_acoustic_rk(:,4))
+
+ potential_acoustic(:) = potential_acoustic_initial_rk(:) + 1.0d0 / 6.0d0 * &
+ (potential_dot_acoustic_rk(:,1) + 2.0d0 * potential_dot_acoustic_rk(:,2) + &
+ 2.0d0 * potential_dot_acoustic_rk(:,3) + potential_dot_acoustic_rk(:,4))
+
+ endif
+
+ endif
+
if(SIMULATION_TYPE ==2)then
b_potential_dot_dot_acoustic = b_potential_dot_dot_acoustic * rmass_inverse_acoustic
b_potential_dot_acoustic = b_potential_dot_acoustic + b_deltatover2*b_potential_dot_dot_acoustic
@@ -4834,7 +4915,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
+ 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
@@ -5441,6 +5522,18 @@
endif
if(time_stepping_scheme == 3)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,:)
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
@@ -5448,18 +5541,6 @@
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,:)
@@ -5472,51 +5553,39 @@
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)
+ 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)
+ 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))
+ (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))
+ (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))
+ (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))
+ (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))
+ (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))
+ (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
More information about the CIG-COMMITS
mailing list