[cig-commits] r19481 - in seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk: EXAMPLES/acoustic_poroelastic EXAMPLES/semi_infinite_homo src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Thu Jan 26 07:45:04 PST 2012
Author: xie.zhinan
Date: 2012-01-26 07:45:03 -0800 (Thu, 26 Jan 2012)
New Revision: 19481
Modified:
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90
Log:
add RK and LDDRK for poro case without consideration of attenuation
Modified: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic 2012-01-26 10:20:11 UTC (rev 19480)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic 2012-01-26 15:45:03 UTC (rev 19481)
@@ -23,10 +23,10 @@
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 5000 # total number of time steps
-deltat = 3.d-4 # duration of a time step
+nt = 3000 # total number of time steps
+deltat = 4.0d-4 # 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 = 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/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D 2012-01-26 10:20:11 UTC (rev 19480)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D 2012-01-26 15:45:03 UTC (rev 19481)
@@ -24,9 +24,9 @@
# time step parameters
nt = 1500 # total number of time steps
-deltat = 2.2e-3 # duration of a time step
+deltat = 1.9e-3 # duration of a time step
USER_T0 = 0.15d0 # 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/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90 2012-01-26 10:20:11 UTC (rev 19480)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90 2012-01-26 15:45:03 UTC (rev 19481)
@@ -448,10 +448,18 @@
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
velocs_poroelastic_LDDRK,displs_poroelastic_LDDRK
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
+ velocs_poroelastic_initial_rk,displs_poroelastic_initial_rk
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
+ accels_poroelastic_rk,velocs_poroelastic_rk
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
accelw_poroelastic,velocw_poroelastic,displw_poroelastic
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
velocw_poroelastic_LDDRK,displw_poroelastic_LDDRK
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
+ velocw_poroelastic_initial_rk,displw_poroelastic_initial_rk
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
+ accelw_poroelastic_rk,velocw_poroelastic_rk
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
accels_poroelastic_adj_coupling, accelw_poroelastic_adj_coupling
double precision, dimension(:), allocatable :: porosity,tortuosity
double precision, dimension(:,:), allocatable :: density,permeability
@@ -2486,6 +2494,17 @@
allocate(velocw_poroelastic_LDDRK(NDIM,nglob_poroelastic))
endif
+ if(time_stepping_scheme == 3)then
+ allocate(accels_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
+ allocate(velocs_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
+ allocate(accelw_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
+ allocate(velocw_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
+ allocate(displs_poroelastic_initial_rk(NDIM,nglob_poroelastic))
+ allocate(velocs_poroelastic_initial_rk(NDIM,nglob_poroelastic))
+ allocate(displw_poroelastic_initial_rk(NDIM,nglob_poroelastic))
+ allocate(velocw_poroelastic_initial_rk(NDIM,nglob_poroelastic))
+ endif
+
! extra array if adjoint and kernels calculation
if(SIMULATION_TYPE == 2 .and. any_poroelastic) then
allocate(b_displs_poroelastic(NDIM,nglob))
@@ -3039,6 +3058,21 @@
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
+
+ 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
+
+ velocw_poroelastic_initial_rk = 0._CUSTOM_REAL
+ displw_poroelastic_initial_rk = 0._CUSTOM_REAL
+
+ endif
+
potential_acoustic = 0._CUSTOM_REAL
potential_dot_acoustic = 0._CUSTOM_REAL
potential_dot_dot_acoustic = 0._CUSTOM_REAL
@@ -6276,7 +6310,6 @@
if(any_poroelastic) then
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(time_stepping_scheme == 1)then
accels_poroelastic(1,:) = accels_poroelastic(1,:) * rmass_s_inverse_poroelastic(:)
@@ -6310,20 +6343,24 @@
displw_poroelastic = displw_poroelastic + beta_LDDRK(i_stage) * displw_poroelastic_LDDRK
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
+ accels_poroelastic(1,:) = accels_poroelastic(1,:) * rmass_s_inverse_poroelastic(:)
+ accels_poroelastic(2,:) = accels_poroelastic(2,:) * rmass_s_inverse_poroelastic(:)
- 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,:)
+ accels_poroelastic_rk(1,:,i_stage) = deltat * accels_poroelastic(1,:)
+ accels_poroelastic_rk(2,:,i_stage) = deltat * accels_poroelastic(2,:)
+ velocs_poroelastic_rk(1,:,i_stage) = deltat * velocs_poroelastic(1,:)
+ velocs_poroelastic_rk(2,:,i_stage) = deltat * velocs_poroelastic(2,:)
- 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,:)
+ accelw_poroelastic(1,:) = accelw_poroelastic(1,:) * rmass_w_inverse_poroelastic(:)
+ accelw_poroelastic(2,:) = accelw_poroelastic(2,:) * rmass_w_inverse_poroelastic(:)
+
+ accelw_poroelastic_rk(1,:,i_stage) = deltat * accelw_poroelastic(1,:)
+ accelw_poroelastic_rk(2,:,i_stage) = deltat * accelw_poroelastic(2,:)
+ velocw_poroelastic_rk(1,:,i_stage) = deltat * velocw_poroelastic(1,:)
+ velocw_poroelastic_rk(2,:,i_stage) = deltat * velocw_poroelastic(2,:)
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
@@ -6333,54 +6370,67 @@
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,:)
+ velocs_poroelastic_initial_rk(1,:) = velocs_poroelastic(1,:)
+ velocs_poroelastic_initial_rk(2,:) = velocs_poroelastic(2,:)
+ displs_poroelastic_initial_rk(1,:) = displs_poroelastic(1,:)
+ displs_poroelastic_initial_rk(2,:) = displs_poroelastic(2,:)
- displ_elastic_initial_rk(1,:) = displ_elastic(1,:)
- displ_elastic_initial_rk(2,:) = displ_elastic(2,:)
- displ_elastic_initial_rk(3,:) = displ_elastic(3,:)
+ velocw_poroelastic_initial_rk(1,:) = velocw_poroelastic(1,:)
+ velocw_poroelastic_initial_rk(2,:) = velocw_poroelastic(2,:)
+ displw_poroelastic_initial_rk(1,:) = displw_poroelastic(1,:)
+ displw_poroelastic_initial_rk(2,:) = displw_poroelastic(2,:)
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)
+ velocs_poroelastic(1,:) = velocs_poroelastic_initial_rk(1,:) + weight_rk * accels_poroelastic_rk(1,:,i_stage)
+ velocs_poroelastic(2,:) = velocs_poroelastic_initial_rk(2,:) + weight_rk * accels_poroelastic_rk(2,:,i_stage)
+ displs_poroelastic(1,:) = displs_poroelastic_initial_rk(1,:) + weight_rk * velocs_poroelastic_rk(1,:,i_stage)
+ displs_poroelastic(2,:) = displs_poroelastic_initial_rk(2,:) + weight_rk * velocs_poroelastic_rk(2,:,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)
+ velocw_poroelastic(1,:) = velocw_poroelastic_initial_rk(1,:) + weight_rk * accelw_poroelastic_rk(1,:,i_stage)
+ velocw_poroelastic(2,:) = velocw_poroelastic_initial_rk(2,:) + weight_rk * accelw_poroelastic_rk(2,:,i_stage)
+ displw_poroelastic(1,:) = displw_poroelastic_initial_rk(1,:) + weight_rk * velocw_poroelastic_rk(1,:,i_stage)
+ displw_poroelastic(2,:) = displw_poroelastic_initial_rk(2,:) + weight_rk * velocw_poroelastic_rk(2,:,i_stage)
+
elseif(i_stage==4)then
- 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))
+ velocs_poroelastic(1,:) = velocs_poroelastic_initial_rk(1,:) + 1.0d0 / 6.0d0 * &
+ (accels_poroelastic_rk(1,:,1) + 2.0d0 * accels_poroelastic_rk(1,:,2) + &
+ 2.0d0 * accels_poroelastic_rk(1,:,3) + accels_poroelastic_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))
+ velocs_poroelastic(2,:) = velocs_poroelastic_initial_rk(2,:) + 1.0d0 / 6.0d0 * &
+ (accels_poroelastic_rk(2,:,1) + 2.0d0 * accels_poroelastic_rk(2,:,2) + &
+ 2.0d0 * accels_poroelastic_rk(2,:,3) + accels_poroelastic_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))
+ displs_poroelastic(1,:) = displs_poroelastic_initial_rk(1,:) + 1.0d0 / 6.0d0 * &
+ (velocs_poroelastic_rk(1,:,1) + 2.0d0 * velocs_poroelastic_rk(1,:,2) + &
+ 2.0d0 * velocs_poroelastic_rk(1,:,3) + velocs_poroelastic_rk(1,:,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))
+ displs_poroelastic(2,:) = displs_poroelastic_initial_rk(2,:) + 1.0d0 / 6.0d0 * &
+ (velocs_poroelastic_rk(2,:,1) + 2.0d0 * velocs_poroelastic_rk(2,:,2) + &
+ 2.0d0 * velocs_poroelastic_rk(2,:,3) + velocs_poroelastic_rk(2,:,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))
+ velocw_poroelastic(1,:) = velocw_poroelastic_initial_rk(1,:) + 1.0d0 / 6.0d0 * &
+ (accelw_poroelastic_rk(1,:,1) + 2.0d0 * accelw_poroelastic_rk(1,:,2) + &
+ 2.0d0 * accelw_poroelastic_rk(1,:,3) + accelw_poroelastic_rk(1,:,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))
+ velocw_poroelastic(2,:) = velocw_poroelastic_initial_rk(2,:) + 1.0d0 / 6.0d0 * &
+ (accelw_poroelastic_rk(2,:,1) + 2.0d0 * accelw_poroelastic_rk(2,:,2) + &
+ 2.0d0 * accelw_poroelastic_rk(2,:,3) + accelw_poroelastic_rk(2,:,4))
+ displw_poroelastic(1,:) = displw_poroelastic_initial_rk(1,:) + 1.0d0 / 6.0d0 * &
+ (velocw_poroelastic_rk(1,:,1) + 2.0d0 * velocw_poroelastic_rk(1,:,2) + &
+ 2.0d0 * velocw_poroelastic_rk(1,:,3) + velocw_poroelastic_rk(1,:,4))
+
+ displw_poroelastic(2,:) = displw_poroelastic_initial_rk(2,:) + 1.0d0 / 6.0d0 * &
+ (velocw_poroelastic_rk(2,:,1) + 2.0d0 * velocw_poroelastic_rk(2,:,2) + &
+ 2.0d0 * velocw_poroelastic_rk(2,:,3) + velocw_poroelastic_rk(2,:,4))
+
endif
endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(SIMULATION_TYPE == 2) then
b_accels_poroelastic(1,:) = b_accels_poroelastic(1,:) * rmass_s_inverse_poroelastic(:)
@@ -6417,6 +6467,9 @@
icount(iglob) = icount(iglob) + 1
if(icount(iglob) ==1)then
+
+ if(time_stepping_scheme == 1)then
+
veloc_elastic(1,iglob) = veloc_elastic(1,iglob) - deltatover2*accel_elastic(1,iglob)
veloc_elastic(3,iglob) = veloc_elastic(3,iglob) - deltatover2*accel_elastic(3,iglob)
accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
@@ -6444,6 +6497,227 @@
velocw_poroelastic(1,iglob) = ZERO
velocw_poroelastic(2,iglob) = ZERO
+ endif
+
+ if(time_stepping_scheme == 2)then
+ ! recovering original velocities and accelerations on boundaries (elastic side)
+ veloc_elastic = veloc_elastic - beta_LDDRK(i_stage) * veloc_elastic_LDDRK
+ displ_elastic = displ_elastic - beta_LDDRK(i_stage) * displ_elastic_LDDRK
+ veloc_elastic_LDDRK = (veloc_elastic_LDDRK - deltat * accel_elastic) / alpha_LDDRK(i_stage)
+ displ_elastic_LDDRK = (displ_elastic_LDDRK - deltat * veloc_elastic) / alpha_LDDRK(i_stage)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
+
+ ! recovering original velocities and accelerations on boundaries (poro side)
+ velocs_poroelastic = velocs_poroelastic - beta_LDDRK(i_stage) * velocs_poroelastic_LDDRK
+ displs_poroelastic = displs_poroelastic - beta_LDDRK(i_stage) * displs_poroelastic_LDDRK
+ velocs_poroelastic_LDDRK = (velocs_poroelastic_LDDRK - deltat * accels_poroelastic) / alpha_LDDRK(i_stage)
+ displs_poroelastic_LDDRK = (velocs_poroelastic_LDDRK - deltat * velocs_poroelastic) / alpha_LDDRK(i_stage)
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) / rmass_s_inverse_poroelastic(iglob)
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
+
+ ! assembling accelerations
+ accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
+ ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ accel_elastic(3,iglob) = ( accel_elastic(3,iglob) + accels_poroelastic(2,iglob) ) / &
+ ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
+ accels_poroelastic(2,iglob) = accel_elastic(3,iglob)
+
+ ! updating velocities
+ ! updating velocities(elastic side)
+ 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
+ ! updating velocities(poro side)
+ velocs_poroelastic_LDDRK = alpha_LDDRK(i_stage) * velocs_poroelastic_LDDRK + deltat * accels_poroelastic
+ displs_poroelastic_LDDRK = alpha_LDDRK(i_stage) * displs_poroelastic_LDDRK + deltat * velocs_poroelastic
+ velocs_poroelastic = velocs_poroelastic + beta_LDDRK(i_stage) * velocs_poroelastic_LDDRK
+ displs_poroelastic = displs_poroelastic + beta_LDDRK(i_stage) * displs_poroelastic_LDDRK
+
+ ! zeros w
+ accelw_poroelastic(1,iglob) = ZERO
+ accelw_poroelastic(2,iglob) = ZERO
+ velocw_poroelastic(1,iglob) = ZERO
+ velocw_poroelastic(2,iglob) = ZERO
+ endif
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if(time_stepping_scheme == 3)then
+
+ ! recovering original velocities and accelerations on boundaries (elastic side)
+ 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
+
+ veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) - weight_rk * accel_elastic_rk(1,iglob,i_stage)
+ veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) - weight_rk * accel_elastic_rk(3,iglob,i_stage)
+ displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) - weight_rk * veloc_elastic_rk(1,iglob,i_stage)
+ displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) - weight_rk * veloc_elastic_rk(3,iglob,i_stage)
+
+
+ elseif(i_stage==4)then
+
+ veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
+ (accel_elastic_rk(1,iglob,1) + 2.0d0 * accel_elastic_rk(1,iglob,2) + &
+ 2.0d0 * accel_elastic_rk(1,iglob,3) + accel_elastic_rk(1,iglob,4))
+
+ veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) - 1.0d0 / 6.0d0 * &
+ (accel_elastic_rk(3,iglob,1) + 2.0d0 * accel_elastic_rk(3,iglob,2) + &
+ 2.0d0 * accel_elastic_rk(3,iglob,3) + accel_elastic_rk(3,iglob,4))
+
+ displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
+ (veloc_elastic_rk(1,iglob,1) + 2.0d0 * veloc_elastic_rk(1,iglob,2) + &
+ 2.0d0 * veloc_elastic_rk(1,iglob,3) + veloc_elastic_rk(1,iglob,4))
+
+ displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) - 1.0d0 / 6.0d0 * &
+ (veloc_elastic_rk(3,iglob,1) + 2.0d0 * veloc_elastic_rk(3,iglob,2) + &
+ 2.0d0 * veloc_elastic_rk(3,iglob,3) + veloc_elastic_rk(3,iglob,4))
+
+ endif
+
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
+
+ accel_elastic_rk(1,iglob,i_stage) = accel_elastic(1,iglob) / deltat
+ accel_elastic_rk(3,iglob,i_stage) = accel_elastic(3,iglob) / deltat
+ veloc_elastic_rk(1,iglob,i_stage) = veloc_elastic(1,iglob) / deltat
+ veloc_elastic_rk(3,iglob,i_stage) = veloc_elastic(3,iglob) / deltat
+
+
+ ! recovering original velocities and accelerations on boundaries (poro side)
+ 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
+
+ velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) - weight_rk * accels_poroelastic_rk(1,iglob,i_stage)
+ velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) - weight_rk * accels_poroelastic_rk(2,iglob,i_stage)
+ displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) - weight_rk * velocs_poroelastic_rk(1,iglob,i_stage)
+ displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) - weight_rk * velocs_poroelastic_rk(2,iglob,i_stage)
+
+
+ elseif(i_stage==4)then
+
+ velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
+ (accels_poroelastic_rk(1,iglob,1) + 2.0d0 * accels_poroelastic_rk(1,iglob,2) + &
+ 2.0d0 * accels_poroelastic_rk(1,iglob,3) + accels_poroelastic_rk(1,iglob,4))
+
+ velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) - 1.0d0 / 6.0d0 * &
+ (accels_poroelastic_rk(2,iglob,1) + 2.0d0 * accels_poroelastic_rk(2,iglob,2) + &
+ 2.0d0 * accels_poroelastic_rk(2,iglob,3) + accels_poroelastic_rk(2,iglob,4))
+
+ displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
+ (velocs_poroelastic_rk(1,iglob,1) + 2.0d0 * velocs_poroelastic_rk(1,iglob,2) + &
+ 2.0d0 * velocs_poroelastic_rk(1,iglob,3) + velocs_poroelastic_rk(1,iglob,4))
+
+ displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) - 1.0d0 / 6.0d0 * &
+ (velocs_poroelastic_rk(2,iglob,1) + 2.0d0 * velocs_poroelastic_rk(2,iglob,2) + &
+ 2.0d0 * velocs_poroelastic_rk(2,iglob,3) + velocs_poroelastic_rk(2,iglob,4))
+
+ endif
+
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) / rmass_s_inverse_poroelastic(iglob)
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
+
+ accels_poroelastic_rk(1,iglob,i_stage) = accels_poroelastic(1,iglob) / deltat
+ accels_poroelastic_rk(2,iglob,i_stage) = accels_poroelastic(2,iglob) / deltat
+ velocs_poroelastic_rk(1,iglob,i_stage) = velocs_poroelastic(1,iglob) / deltat
+ velocs_poroelastic_rk(2,iglob,i_stage) = velocs_poroelastic(2,iglob) / deltat
+
+
+ ! assembling accelerations
+ accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
+ ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ accel_elastic(3,iglob) = ( accel_elastic(3,iglob) + accels_poroelastic(2,iglob) ) / &
+ ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
+ accels_poroelastic(2,iglob) = accel_elastic(3,iglob)
+
+ ! updating velocities
+ ! updating velocities(elastic side)
+
+ accel_elastic_rk(1,iglob,i_stage) = accel_elastic(1,iglob) * deltat
+ accel_elastic_rk(3,iglob,i_stage) = accel_elastic(3,iglob) * deltat
+
+ 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
+
+ veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) + weight_rk * accel_elastic_rk(1,iglob,i_stage)
+ veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) + weight_rk * accel_elastic_rk(3,iglob,i_stage)
+ displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) + weight_rk * veloc_elastic_rk(1,iglob,i_stage)
+ displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) + weight_rk * veloc_elastic_rk(3,iglob,i_stage)
+
+
+ elseif(i_stage==4)then
+
+ veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
+ (accel_elastic_rk(1,iglob,1) + 2.0d0 * accel_elastic_rk(1,iglob,2) + &
+ 2.0d0 * accel_elastic_rk(1,iglob,3) + accel_elastic_rk(1,iglob,4))
+
+ veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) + 1.0d0 / 6.0d0 * &
+ (accel_elastic_rk(3,iglob,1) + 2.0d0 * accel_elastic_rk(3,iglob,2) + &
+ 2.0d0 * accel_elastic_rk(3,iglob,3) + accel_elastic_rk(3,iglob,4))
+
+ displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
+ (veloc_elastic_rk(1,iglob,1) + 2.0d0 * veloc_elastic_rk(1,iglob,2) + &
+ 2.0d0 * veloc_elastic_rk(1,iglob,3) + veloc_elastic_rk(1,iglob,4))
+
+ displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) + 1.0d0 / 6.0d0 * &
+ (veloc_elastic_rk(3,iglob,1) + 2.0d0 * veloc_elastic_rk(3,iglob,2) + &
+ 2.0d0 * veloc_elastic_rk(3,iglob,3) + veloc_elastic_rk(3,iglob,4))
+
+ endif
+ ! updating velocities(poro side)
+
+ accels_poroelastic_rk(1,iglob,i_stage) = deltat * accels_poroelastic(1,iglob)
+ accels_poroelastic_rk(2,iglob,i_stage) = deltat * accels_poroelastic(2,iglob)
+ velocs_poroelastic_rk(1,iglob,i_stage) = deltat * velocs_poroelastic(1,iglob)
+ velocs_poroelastic_rk(2,iglob,i_stage) = deltat * velocs_poroelastic(2,iglob)
+
+
+ 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
+
+ velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) + weight_rk * accels_poroelastic_rk(1,iglob,i_stage)
+ velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) + weight_rk * accels_poroelastic_rk(2,iglob,i_stage)
+ displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) + weight_rk * velocs_poroelastic_rk(1,iglob,i_stage)
+ displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) + weight_rk * velocs_poroelastic_rk(2,iglob,i_stage)
+
+
+ elseif(i_stage==4)then
+
+ velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
+ (accels_poroelastic_rk(1,iglob,1) + 2.0d0 * accels_poroelastic_rk(1,iglob,2) + &
+ 2.0d0 * accels_poroelastic_rk(1,iglob,3) + accels_poroelastic_rk(1,iglob,4))
+
+ velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) + 1.0d0 / 6.0d0 * &
+ (accels_poroelastic_rk(2,iglob,1) + 2.0d0 * accels_poroelastic_rk(2,iglob,2) + &
+ 2.0d0 * accels_poroelastic_rk(2,iglob,3) + accels_poroelastic_rk(2,iglob,4))
+
+ displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
+ (velocs_poroelastic_rk(1,iglob,1) + 2.0d0 * velocs_poroelastic_rk(1,iglob,2) + &
+ 2.0d0 * velocs_poroelastic_rk(1,iglob,3) + velocs_poroelastic_rk(1,iglob,4))
+
+ displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) + 1.0d0 / 6.0d0 * &
+ (velocs_poroelastic_rk(2,iglob,1) + 2.0d0 * velocs_poroelastic_rk(2,iglob,2) + &
+ 2.0d0 * velocs_poroelastic_rk(2,iglob,3) + velocs_poroelastic_rk(2,iglob,4))
+
+ endif
+
+ endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
if(SIMULATION_TYPE == 2) then
b_veloc_elastic(1,iglob) = b_veloc_elastic(1,iglob) - b_deltatover2*b_accel_elastic(1,iglob)
b_veloc_elastic(3,iglob) = b_veloc_elastic(3,iglob) - b_deltatover2*b_accel_elastic(3,iglob)
More information about the CIG-COMMITS
mailing list