[cig-commits] r19426 - 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 04:02:12 PST 2012
Author: xie.zhinan
Date: 2012-01-24 04:02:12 -0800 (Tue, 24 Jan 2012)
New Revision: 19426
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_poro_fluid.f90
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_poro_solid.f90
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/prepare_source_time_function.f90
seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90
Log:
change source time function for implementing other time scheme
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 08:54:58 UTC (rev 19425)
+++ 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)
@@ -16,7 +16,7 @@
add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
-ATTENUATION_VISCOELASTIC_SOLID = .true. # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_VISCOELASTIC_SOLID = .false. # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
ATTENUATION_PORO_FLUID_PART = .false. # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
Q0 = 1 # quality factor for viscous attenuation
freq0 = 10 # frequency for viscous attenuation
@@ -26,7 +26,7 @@
nt = 1500 # total number of time steps
deltat = 7.5e-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 = 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 = 1 # 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_poro_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_poro_fluid.f90 2012-01-24 08:54:58 UTC (rev 19425)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_poro_fluid.f90 2012-01-24 12:02:12 UTC (rev 19426)
@@ -60,7 +60,8 @@
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
C_k,M_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
b_absorb_poro_w_left,b_absorb_poro_w_right,b_absorb_poro_w_bottom,b_absorb_poro_w_top,&
- nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,f0,freq0,Q0)
+ nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,f0,freq0,Q0, &
+ stage_time_scheme,i_stage) !xiezhinan
! compute forces for the fluid poroelastic part
@@ -77,6 +78,7 @@
integer, dimension(nelemabs) :: ib_right
integer, dimension(nelemabs) :: ib_bottom
integer, dimension(nelemabs) :: ib_top
+ integer :: stage_time_scheme,i_stage !xiezhinan
logical :: anyabs,initialfield,ATTENUATION_VISCOELASTIC_SOLID
logical :: SAVE_FORWARD
@@ -99,7 +101,7 @@
double precision, dimension(numat) :: porosity,tortuosity
double precision, dimension(4,3,numat) :: poroelastcoef
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
- real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP) :: source_time_function
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP,stage_time_scheme) :: source_time_function !xiezhinan
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(nglob) :: C_k,M_k
@@ -832,7 +834,7 @@
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_source(i_source))
accelw_poroelastic(:,iglob) = accelw_poroelastic(:,iglob) + &
- (1._CUSTOM_REAL - rhol_f/rhol_bar)*sourcearray(i_source,:,i,j)*source_time_function(i_source,it)
+ (1._CUSTOM_REAL - rhol_f/rhol_bar)*sourcearray(i_source,:,i,j)*source_time_function(i_source,it,i_stage)
enddo
enddo
else ! backward wavefield
@@ -840,7 +842,8 @@
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_source(i_source))
b_accelw_poroelastic(:,iglob) = b_accelw_poroelastic(:,iglob) + &
- (1._CUSTOM_REAL - rhol_f/rhol_bar)*sourcearray(i_source,:,i,j)*source_time_function(i_source,NSTEP-it+1)
+ (1._CUSTOM_REAL - rhol_f/rhol_bar)*sourcearray(i_source,:,i,j) * &
+ source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
enddo
enddo
endif !endif SIMULATION_TYPE == 1
Modified: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_poro_solid.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_poro_solid.f90 2012-01-24 08:54:58 UTC (rev 19425)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_poro_solid.f90 2012-01-24 12:02:12 UTC (rev 19426)
@@ -60,7 +60,8 @@
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
mufr_k,B_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
b_absorb_poro_s_left,b_absorb_poro_s_right,b_absorb_poro_s_bottom,b_absorb_poro_s_top,&
- nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,f0,freq0,Q0)
+ nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,f0,freq0,Q0,&
+ stage_time_scheme,i_stage) !xiezhinan
! compute forces for the solid poroelastic part
@@ -77,6 +78,7 @@
integer, dimension(nelemabs) :: ib_right
integer, dimension(nelemabs) :: ib_bottom
integer, dimension(nelemabs) :: ib_top
+ integer :: stage_time_scheme,i_stage !xiezhinan
logical :: anyabs,initialfield,ATTENUATION_VISCOELASTIC_SOLID
logical :: SAVE_FORWARD
@@ -99,7 +101,7 @@
double precision, dimension(numat) :: porosity,tortuosity
double precision, dimension(4,3,numat) :: poroelastcoef
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
- real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP) :: source_time_function
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP,stage_time_scheme) :: source_time_function !xiezhinan
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(nglob) :: mufr_k,B_k
@@ -839,7 +841,7 @@
do i_source=1,NSOURCES
! if this processor core carries the source and the source element is poroelastic
- if (is_proc_source(i_source) == 1 .and. poroelastic(ispec_selected_source(i_source))) then
+ if (is_proc_source(i_source) == 1 .and. poroelastic(ispec_selected_source(i_source))) then
phil = porosity(kmato(ispec_selected_source(i_source)))
tortl = tortuosity(kmato(ispec_selected_source(i_source)))
@@ -853,7 +855,7 @@
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_source(i_source))
accels_poroelastic(:,iglob) = accels_poroelastic(:,iglob) + &
- (1._CUSTOM_REAL - phil/tortl)*sourcearray(i_source,:,i,j)*source_time_function(i_source,it)
+ (1._CUSTOM_REAL - phil/tortl)*sourcearray(i_source,:,i,j)*source_time_function(i_source,it,i_stage)
enddo
enddo
else ! backward wavefield
@@ -861,7 +863,8 @@
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_source(i_source))
b_accels_poroelastic(:,iglob) = b_accels_poroelastic(:,iglob) + &
- (1._CUSTOM_REAL - phil/tortl)*sourcearray(i_source,:,i,j)*source_time_function(i_source,NSTEP-it+1)
+ (1._CUSTOM_REAL - phil/tortl)*sourcearray(i_source,:,i,j) * &
+ source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
enddo
enddo
endif !endif SIMULATION_TYPE == 1
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-24 08:54:58 UTC (rev 19425)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2012-01-24 12:02:12 UTC (rev 19426)
@@ -59,7 +59,8 @@
v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot,&
nleft,nright,nbot,over_critical_angle,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,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)
+ nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k,&
+ stage_time_scheme,i_stage) !xiezhinan
! compute forces for the elastic elements
@@ -79,6 +80,7 @@
integer, dimension(nelemabs) :: ib_right
integer, dimension(nelemabs) :: ib_bottom
integer, dimension(nelemabs) :: ib_top
+ integer :: stage_time_scheme,i_stage !xiezhinan
logical :: anyabs,assign_external_model,initialfield,ATTENUATION_VISCOELASTIC_SOLID,add_Bielak_conditions
@@ -102,7 +104,7 @@
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,vsext,rhoext
double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
- real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP) :: source_time_function
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP,stage_time_scheme) :: source_time_function !xiezhinan
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_displ_elastic
@@ -899,9 +901,9 @@
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_source(i_source))
accel_elastic(1,iglob) = accel_elastic(1,iglob) + &
- sourcearray(i_source,1,i,j)*source_time_function(i_source,it)
+ sourcearray(i_source,1,i,j)*source_time_function(i_source,it,i_stage)
accel_elastic(3,iglob) = accel_elastic(3,iglob) + &
- sourcearray(i_source,2,i,j)*source_time_function(i_source,it)
+ sourcearray(i_source,2,i,j)*source_time_function(i_source,it,i_stage)
enddo
enddo
else ! backward wavefield
@@ -909,9 +911,9 @@
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_source(i_source))
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + &
- sourcearray(i_source,1,i,j)*source_time_function(i_source,NSTEP-it+1)
+ sourcearray(i_source,1,i,j)*source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) + &
- sourcearray(i_source,2,i,j)*source_time_function(i_source,NSTEP-it+1)
+ sourcearray(i_source,2,i,j)*source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
enddo
enddo
endif !endif SIMULATION_TYPE == 1
Modified: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/prepare_source_time_function.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/prepare_source_time_function.f90 2012-01-24 08:54:58 UTC (rev 19425)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/prepare_source_time_function.f90 2012-01-24 12:02:12 UTC (rev 19426)
@@ -46,7 +46,7 @@
subroutine prepare_source_time_function(myrank,NSTEP,NSOURCES,source_time_function, &
time_function_type,f0,tshift_src,factor,aval, &
- t0,nb_proc_source,deltat)
+ t0,nb_proc_source,deltat,stage_time_scheme,c_LDDRK) !xiezhinan
! prepares source_time_function array
@@ -62,15 +62,25 @@
double precision :: t0
integer,dimension(NSOURCES) :: nb_proc_source
double precision :: deltat
+ integer :: stage_time_scheme !xiezhinan
+ double precision, dimension(stage) :: c_LDDRK !xiezhinan
- real(kind=CUSTOM_REAL),dimension(NSOURCES,NSTEP) :: source_time_function
+ real(kind=CUSTOM_REAL),dimension(NSOURCES,NSTEP,stage_time_scheme) :: source_time_function
! local parameters
double precision :: stf_used,time
double precision, dimension(NSOURCES) :: hdur,hdur_gauss
double precision, external :: netlib_specfun_erf
integer :: it,i_source
+ integer :: i_stage
+ double precision, dimension(4) :: c_RK
+ if(stage_time_scheme == 4)then
+ c_RK(1)=0.0d0*deltat
+ c_RK(2)=0.5d0*deltat
+ c_RK(3)=0.5d0*deltat
+ c_RK(4)=1.0d0*deltat
+ endif
! user output
if (myrank == 0) then
@@ -84,14 +94,28 @@
! do i_source=1,NSOURCES
! loop on all the time steps
+
+
do it = 1,NSTEP
! note: t0 is the simulation start time, tshift_src is the time shift of the source
! relative to this start time
+ do i_stage = 1,stage_time_scheme
! compute current time
- time = (it-1)*deltat
+ if(stage_time_scheme == 1)then
+ time = (it-1)*deltat
+ endif
+ if(stage_time_scheme == 4)then
+ time = (it-1)*deltat+c_RK(i_stage)*deltat
+ endif
+
+ if(stage_time_scheme == 6)then
+ time = (it-1)*deltat+c_LDDRK(i_stage)*deltat
+ endif
+
+
stf_used = 0.d0
! loop on all the sources
@@ -100,7 +124,7 @@
if( time_function_type(i_source) == 1 ) then
! Ricker (second derivative of a Gaussian) source time function
- source_time_function(i_source,it) = - factor(i_source) * &
+ source_time_function(i_source,it,i_stage) = - factor(i_source) * &
(ONE-TWO*aval(i_source)*(time-t0-tshift_src(i_source))**2) * &
exp(-aval(i_source)*(time-t0-tshift_src(i_source))**2)
@@ -111,14 +135,14 @@
else if( time_function_type(i_source) == 2 ) then
! first derivative of a Gaussian source time function
- source_time_function(i_source,it) = - factor(i_source) * &
+ source_time_function(i_source,it,i_stage) = - factor(i_source) * &
TWO*aval(i_source)*(time-t0-tshift_src(i_source)) * &
exp(-aval(i_source)*(time-t0-tshift_src(i_source))**2)
else if(time_function_type(i_source) == 3 .or. time_function_type(i_source) == 4) then
! Gaussian or Dirac (we use a very thin Gaussian instead) source time function
- source_time_function(i_source,it) = factor(i_source) * &
+ source_time_function(i_source,it,i_stage) = factor(i_source) * &
exp(-aval(i_source)*(time-t0-tshift_src(i_source))**2)
else if(time_function_type(i_source) == 5) then
@@ -126,14 +150,14 @@
! Heaviside source time function (we use a very thin error function instead)
hdur(i_source) = 1.d0 / f0(i_source)
hdur_gauss(i_source) = hdur(i_source) * 5.d0 / 3.d0
- source_time_function(i_source,it) = factor(i_source) * 0.5d0*(1.0d0 + &
+ source_time_function(i_source,it,i_stage) = factor(i_source) * 0.5d0*(1.0d0 + &
netlib_specfun_erf(SOURCE_DECAY_MIMIC_TRIANGLE*(time-t0-tshift_src(i_source))/hdur_gauss(i_source)))
else
call exit_MPI('unknown source time function')
endif
- stf_used = stf_used + source_time_function(i_source,it)
+ stf_used = stf_used + source_time_function(i_source,it,i_stage)
enddo
@@ -144,7 +168,8 @@
write(55,*) sngl(time-t0),sngl(stf_used),sngl(time)
endif
- !enddo
+ enddo
+
enddo
if (myrank == 0) close(55)
@@ -154,7 +179,7 @@
! since source contribution is linear, the source_time_function is cut down by that number (it would have been similar
! if we just had elected one of those processes).
do i_source=1,NSOURCES
- source_time_function(i_source,:) = source_time_function(i_source,:) / nb_proc_source(i_source)
+ source_time_function(i_source,:,:) = source_time_function(i_source,:,:) / nb_proc_source(i_source)
enddo
end subroutine prepare_source_time_function
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 08:54:58 UTC (rev 19425)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/src/specfem2D/specfem2D.F90 2012-01-24 12:02:12 UTC (rev 19426)
@@ -435,6 +435,7 @@
double precision :: mul_unrelaxed_elastic,lambdal_unrelaxed_elastic,lambdaplus2mu_unrelaxed_elastic,kappal
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_adj_coupling
double precision, dimension(:,:), allocatable :: &
coord, flagrange,xinterp,zinterp,Uxinterp,Uzinterp,vector_field_display
@@ -494,7 +495,7 @@
integer, dimension(:), allocatable :: ispec_selected_source,iglob_source,&
is_proc_source,nb_proc_source
double precision, dimension(:), allocatable :: aval
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: source_time_function
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: source_time_function
double precision, external :: netlib_specfun_erf
double precision :: vpImin,vpImax,vpIImin,vpIImax
@@ -889,9 +890,7 @@
logical :: save_everywhere = .false.
! for LDDRK46
- logical :: LDDRK
- integer :: i_LDDRK
- double precision :: time_LDDRK
+ integer :: i_stage,stage_time_scheme
double precision, dimension(6):: alpha_LDDRK,beta_LDDRK,c_LDDRK
!***********************************************************************
@@ -1001,11 +1000,16 @@
t0,initialfield,ipass,deltat,USER_T0)
- !
- !---- read Par_LDDRK
- !
- if(time_stepping_scheme == 2)then !xiezhinan
- call read_databases_LDDRK(alpha_LDDRK,beta_LDDRK,c_LDDRK) !xiezhinan
+
+ !---- read Par_LDDRK
+
+ if(time_stepping_scheme == 1)then
+ stage_time_scheme=1
+ elseif(time_stepping_scheme == 2)then !xiezhinan
+ stage_time_scheme=stage
+ call read_databases_LDDRK(alpha_LDDRK,beta_LDDRK,c_LDDRK) !xiezhinan
+ elseif(time_stepping_scheme == 3)then
+ stage_time_scheme=4
endif !xiezhinan
@@ -2352,6 +2356,11 @@
allocate(accel_elastic(3,nglob_elastic))
allocate(accel_elastic_adj_coupling(3,nglob_elastic))
allocate(rmass_inverse_elastic(nglob_elastic))
+
+ if(time_stepping_scheme==2)then
+ allocate(displ_elastic_LDDRK(3,nglob_elastic))
+ allocate(veloc_elastic_LDDRK(3,nglob_elastic))
+ endif
! extra array if adjoint and kernels calculation
if(SIMULATION_TYPE == 2 .and. any_elastic) then
@@ -2928,6 +2937,11 @@
veloc_elastic = 0._CUSTOM_REAL
accel_elastic = 0._CUSTOM_REAL
+ if(time_stepping_scheme == 2) then !xiezhinan
+ displ_elastic_LDDRK = 0._CUSTOM_REAL
+ veloc_elastic_LDDRK = 0._CUSTOM_REAL
+ endif
+
displs_poroelastic = 0._CUSTOM_REAL
velocs_poroelastic = 0._CUSTOM_REAL
accels_poroelastic = 0._CUSTOM_REAL
@@ -3193,17 +3207,17 @@
! compute the source time function and store it in a text file
if(.not. initialfield) then
- allocate(source_time_function(NSOURCES,NSTEP))
- source_time_function(:,:) = 0._CUSTOM_REAL
+ allocate(source_time_function(NSOURCES,NSTEP,stage_time_scheme)) !xiezhinan
+ source_time_function(:,:,:) = 0._CUSTOM_REAL
! computes source time function array
call prepare_source_time_function(myrank,NSTEP,NSOURCES,source_time_function, &
time_function_type,f0,tshift_src,factor,aval, &
- t0,nb_proc_source,deltat)
+ t0,nb_proc_source,deltat,stage_time_scheme,c_LDDRK) !xiezhinan
else
! uses an initialfield
! dummy allocation
- allocate(source_time_function(1,1))
+ allocate(source_time_function(1,1,1)) !xiezhinan
endif
! determine if coupled fluid-solid simulation
@@ -4144,14 +4158,24 @@
! compute current time
time = (it-1)*deltat
+ do i_stage=1, stage_time_scheme !xiezhinan
+
! 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 &
+ deltatsquareover2*accel_elastic
veloc_elastic = veloc_elastic + deltatover2*accel_elastic
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
if(SIMULATION_TYPE == 2) then ! Adjoint calculation
b_displ_elastic = b_displ_elastic &
@@ -4163,6 +4187,8 @@
endif
if(any_poroelastic) then
+
+ if(time_stepping_scheme==1)then
!for the solid
displs_poroelastic = displs_poroelastic &
+ deltat*velocs_poroelastic &
@@ -4177,7 +4203,22 @@
velocw_poroelastic = velocw_poroelastic + deltatover2*accelw_poroelastic
accelw_poroelastic_adj_coupling = accelw_poroelastic
accelw_poroelastic = ZERO
+ endif
+ if(time_stepping_scheme==2)then
+ !for the solid
+ accels_poroelastic = ZERO
+ !for the fluid
+ accelw_poroelastic = ZERO
+ endif
+
+ if(time_stepping_scheme==3)then
+ !for the solid
+ accels_poroelastic = ZERO
+ !for the fluid
+ accelw_poroelastic = ZERO
+ endif
+
if(SIMULATION_TYPE == 2) then ! Adjoint calculation
!for the solid
b_displs_poroelastic = b_displs_poroelastic &
@@ -4259,7 +4300,7 @@
!-----------------------------------------
if(any_acoustic) then
-
+ if(time_stepping_scheme==1)then
! Newmark time scheme
potential_acoustic = potential_acoustic &
+ deltat*potential_dot_acoustic &
@@ -4267,7 +4308,16 @@
potential_dot_acoustic = potential_dot_acoustic &
+ deltatover2*potential_dot_dot_acoustic
potential_dot_dot_acoustic = ZERO
+ endif
+ if(time_stepping_scheme==2)then
+ potential_dot_dot_acoustic = ZERO
+ endif
+
+ if(time_stepping_scheme==3)then
+ potential_dot_dot_acoustic = ZERO
+ endif
+
if(SIMULATION_TYPE == 2) then ! Adjoint calculation
b_potential_acoustic = b_potential_acoustic &
+ b_deltat*b_potential_dot_acoustic &
@@ -4588,7 +4638,7 @@
iglob = ibool(i,j,ispec_selected_source(i_source))
hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - source_time_function(i_source,it)*hlagrange
+ - source_time_function(i_source,it,i_stage)*hlagrange
enddo
enddo
else
@@ -4598,7 +4648,7 @@
iglob = ibool(i,j,ispec_selected_source(i_source))
hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
- - source_time_function(i_source,NSTEP-it+1)*hlagrange
+ - source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)*hlagrange
enddo
enddo
endif
@@ -4733,7 +4783,8 @@
count_left,count_right,count_bottom,over_critical_angle, &
NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
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)
+ nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k, &
+ stage_time_scheme,i_stage) !xiezhinan)
if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
!--- left absorbing boundary
@@ -5217,9 +5268,9 @@
iglob = ibool(i,j,ispec_selected_source(i_source))
hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
accel_elastic(1,iglob) = accel_elastic(1,iglob) &
- - sin(angleforce(i_source))*source_time_function(i_source,it)*hlagrange
+ - sin(angleforce(i_source))*source_time_function(i_source,it,i_stage)*hlagrange
accel_elastic(3,iglob) = accel_elastic(3,iglob) &
- + cos(angleforce(i_source))*source_time_function(i_source,it)*hlagrange
+ + cos(angleforce(i_source))*source_time_function(i_source,it,i_stage)*hlagrange
enddo
enddo
else ! SH (membrane) calculation
@@ -5228,7 +5279,7 @@
iglob = ibool(i,j,ispec_selected_source(i_source))
hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
accel_elastic(2,iglob) = accel_elastic(2,iglob) &
- + source_time_function(i_source,it)*hlagrange
+ + source_time_function(i_source,it,i_stage)*hlagrange
enddo
enddo
endif
@@ -5239,10 +5290,10 @@
iglob = ibool(i,j,ispec_selected_source(i_source))
hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) &
- - sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1) &
+ - sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1) &
*hlagrange
b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) &
- + cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1) &
+ + cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1) &
*hlagrange
enddo
enddo
@@ -5252,7 +5303,7 @@
iglob = ibool(i,j,ispec_selected_source(i_source))
hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) &
- + source_time_function(i_source,NSTEP-it+1)*hlagrange
+ + source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)*hlagrange
enddo
enddo
endif
@@ -5317,12 +5368,32 @@
! ************************************************************************************
if(any_elastic) 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
- veloc_elastic = veloc_elastic + deltatover2*accel_elastic
+ if(time_stepping_scheme == 1)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
+ veloc_elastic = veloc_elastic + deltatover2*accel_elastic
+ endif
+ if(time_stepping_scheme == 2)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
+
+ 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
+ endif
+
+ if(time_stepping_scheme == 3)then
+
+ 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(:)
@@ -5367,7 +5438,8 @@
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
mufr_k,B_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
b_absorb_poro_s_left,b_absorb_poro_s_right,b_absorb_poro_s_bottom,b_absorb_poro_s_top,&
- nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,f0(1),freq0,Q0)
+ nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,f0(1),freq0,Q0,&
+ stage_time_scheme,i_stage) !xiezhinan)
@@ -5390,7 +5462,8 @@
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
C_k,M_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
b_absorb_poro_w_left,b_absorb_poro_w_right,b_absorb_poro_w_bottom,b_absorb_poro_w_top,&
- nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,f0(1),freq0,Q0)
+ nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,f0(1),freq0,Q0,&
+ stage_time_scheme,i_stage) !xiezhinan)
if(SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
@@ -5901,14 +5974,14 @@
hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
! s
accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) - hlagrange * &
- (1._CUSTOM_REAL - phil/tortl)*sin(angleforce(i_source))*source_time_function(i_source,it)
+ (1._CUSTOM_REAL - phil/tortl)*sin(angleforce(i_source))*source_time_function(i_source,it,i_stage)
accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + hlagrange * &
- (1._CUSTOM_REAL - phil/tortl)*cos(angleforce(i_source))*source_time_function(i_source,it)
+ (1._CUSTOM_REAL - phil/tortl)*cos(angleforce(i_source))*source_time_function(i_source,it,i_stage)
! w
accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) - hlagrange * &
- (1._CUSTOM_REAL - rhol_f/rhol_bar)*sin(angleforce(i_source))*source_time_function(i_source,it)
+ (1._CUSTOM_REAL - rhol_f/rhol_bar)*sin(angleforce(i_source))*source_time_function(i_source,it,i_stage)
accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) + hlagrange * &
- (1._CUSTOM_REAL - rhol_f/rhol_bar)*cos(angleforce(i_source))*source_time_function(i_source,it)
+ (1._CUSTOM_REAL - rhol_f/rhol_bar)*cos(angleforce(i_source))*source_time_function(i_source,it,i_stage)
enddo
enddo
else ! backward wavefield
@@ -5918,14 +5991,18 @@
hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
! b_s
b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) - hlagrange * &
- (1._CUSTOM_REAL - phil/tortl)*sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ (1._CUSTOM_REAL - phil/tortl)*sin(angleforce(i_source))* &
+ source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + hlagrange * &
- (1._CUSTOM_REAL - phil/tortl)*cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ (1._CUSTOM_REAL - phil/tortl)*cos(angleforce(i_source))* &
+ source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
!b_w
b_accelw_poroelastic(1,iglob) = b_accelw_poroelastic(1,iglob) - hlagrange * &
- (1._CUSTOM_REAL - rhol_f/rhol_bar)*sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ (1._CUSTOM_REAL - rhol_f/rhol_bar)*sin(angleforce(i_source))* &
+ source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
b_accelw_poroelastic(2,iglob) = b_accelw_poroelastic(2,iglob) + hlagrange * &
- (1._CUSTOM_REAL - rhol_f/rhol_bar)*cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ (1._CUSTOM_REAL - rhol_f/rhol_bar)*cos(angleforce(i_source))* &
+ source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
enddo
enddo
endif !endif SIMULATION_TYPE == 1
@@ -6072,6 +6149,8 @@
enddo
endif
+ enddo !xiezhinan LDDRK or RK
+
! ********************************************************************************************
! reading lastframe for adjoint/kernels calculation
! ********************************************************************************************
More information about the CIG-COMMITS
mailing list