[cig-commits] r19530 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES/acoustic_poroelastic EXAMPLES/attenuation EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom EXAMPLES/semi_infinite_homo src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Tue Jan 31 04:45:00 PST 2012
Author: xie.zhinan
Date: 2012-01-31 04:44:59 -0800 (Tue, 31 Jan 2012)
New Revision: 19530
Modified:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic
seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/interfaces_acoustic_poroelastic.dat
seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D
seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom/Par_file_fluid_solid
seismo/2D/SPECFEM2D/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
RK and LDDRK time schemes not supported for
adjoint inversion, mixed elastic/poroelastic forword simulations
and poroelastic forword simulations with attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic 2012-01-31 04:55:42 UTC (rev 19529)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/Par_file_acoustic_poroelastic 2012-01-31 12:44:59 UTC (rev 19530)
@@ -19,12 +19,12 @@
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 = 1000 # frequency for viscous attenuation
+freq0 = 1000 # frequency for viscous attenuation
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 = 1.2d-4 # duration of a time step
+nt = 2500 # total number of time steps
+deltat = 4.d-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
@@ -91,8 +91,8 @@
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
-1 3 2650.d0 880.d0 0.3d0 1.2 1d-9 0.d0 1d-9 1.22d10 1.985d9 9.6d9 10d-4 0.0001 10000
-2 1 2500.d0 3400.d0 1963.d0 0 0 9999 9999 0 0 0 0 0 0
+1 3 2650.d0 880.d0 0.3d0 1.2 1d-9 0.d0 1d-9 1.22d10 1.985d9 9.6d9 10d-4 0.0001 9999999999
+2 1 1020.d0 1500.d0 0.d0 0 0 10.d0 10.d0 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -130,7 +130,7 @@
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
xmax = 4800.d0 # abscissa of right side of the model
-nx = 260 # number of elements along X
+nx = 130 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
@@ -140,5 +140,5 @@
# define the different regions of the model in the (nx,nz) spectral element mesh
nbregions = 2 # nb of regions and model number for each
-1 260 1 110 1
-1 260 111 220 2
+1 130 1 55 1
+1 130 56 110 2
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/interfaces_acoustic_poroelastic.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/interfaces_acoustic_poroelastic.dat 2012-01-31 04:55:42 UTC (rev 19529)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/acoustic_poroelastic/interfaces_acoustic_poroelastic.dat 2012-01-31 12:44:59 UTC (rev 19530)
@@ -29,6 +29,6 @@
#
# layer number 1 (bottom layer)
#
- 110
+ 55
#
- 110
+ 55
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D 2012-01-31 04:55:42 UTC (rev 19529)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/Par_file_attenuation_2D 2012-01-31 12:44:59 UTC (rev 19530)
@@ -23,10 +23,10 @@
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 3000 # total number of time steps
-deltat = 5e-4 # duration of a time step
+nt = 1400 # total number of time steps
+deltat = 1e-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/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom/Par_file_fluid_solid
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom/Par_file_fluid_solid 2012-01-31 04:55:42 UTC (rev 19529)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/from_2000_Geophysics_paper_sinusoidal_ocean_bottom/Par_file_fluid_solid 2012-01-31 12:44:59 UTC (rev 19530)
@@ -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 = 0.55d-3 # duration of a time step
+nt = 2000 # total number of time steps
+deltat = 0.7d-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/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D 2012-01-31 04:55:42 UTC (rev 19529)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/semi_infinite_homo/Par_file_elastic_2D 2012-01-31 12:44:59 UTC (rev 19530)
@@ -24,9 +24,9 @@
# time step parameters
nt = 1500 # total number of time steps
-deltat = 1.9e-3 # duration of a time step
+deltat = 0.8e-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 = 3 # 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/trunk/src/specfem2D/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2012-01-31 04:55:42 UTC (rev 19529)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2012-01-31 12:44:59 UTC (rev 19530)
@@ -205,7 +205,7 @@
! compute Grad(displ_elastic) at time step n for attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
- temp_displ_elastic = displ_elastic + deltat * veloc_elastic + 0.5 * deltat**2 * accel_elastic !xiezhinan
+ temp_displ_elastic = displ_elastic + deltat * veloc_elastic + 0.5 * deltat**2 * accel_elastic
call compute_gradient_attenuation(displ_elastic,dux_dxl_n,duz_dxl_n, &
dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
endif
@@ -990,18 +990,18 @@
! evolution e1
if(stage_time_scheme == 1) then
- Un = e1(i,j,ispec,i_sls)
- tauinv = - inv_tau_sigma_nu1(i,j,ispec,i_sls)
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = theta_n * phi_nu1(i,j,ispec,i_sls)
- Snp1 = theta_np1 * phi_nu1(i,j,ispec,i_sls)
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e1(i,j,ispec,i_sls) = Unp1
+ Un = e1(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu1(i,j,ispec,i_sls)
+ tauinvsquare = tauinv * tauinv
+ tauinvcube = tauinvsquare * tauinv
+ tauinvUn = tauinv * Un
+ Sn = theta_n * phi_nu1(i,j,ispec,i_sls)
+ Snp1 = theta_np1 * phi_nu1(i,j,ispec,i_sls)
+ Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e1(i,j,ispec,i_sls) = Unp1
endif
if(stage_time_scheme == 6) then
@@ -1009,158 +1009,129 @@
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_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e1_LDDRK(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_time_scheme - Un * tauinv)
- e1(i,j,ispec,i_sls) = Un + beta_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls)
+ 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)
+ 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 .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
+ if(i_stage==1)then
+ e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
+ endif
- e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
+ 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)
- endif
+ elseif(i_stage==4)then
- 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)
+ 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
- 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))
-
+ ! evolution e11
+ if(stage_time_scheme == 1) then
+ Un = e11(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
+ tauinvsquare = tauinv * tauinv
+ tauinvcube = tauinvsquare * tauinv
+ tauinvUn = tauinv * Un
+ Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i,j,ispec,i_sls)
+ Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i,j,ispec,i_sls)
+ Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e11(i,j,ispec,i_sls) = Unp1
endif
+ if(stage_time_scheme == 6) then
+ 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_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
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! evolution e11
- if(stage_time_scheme == 1) then
- Un = e11(i,j,ispec,i_sls)
- tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i,j,ispec,i_sls)
- Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i,j,ispec,i_sls)
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e11(i,j,ispec,i_sls) = Unp1
- endif
-
- if(stage_time_scheme == 6) then
- 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_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))
+ 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
-
+ 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)
- tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
- Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e13(i,j,ispec,i_sls) = Unp1
- endif
+ if(stage_time_scheme == 1) then
+ Un = e13(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
+ tauinvsquare = tauinv * tauinv
+ tauinvcube = tauinvsquare * tauinv
+ tauinvUn = tauinv * Un
+ Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
+ Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
+ Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e13(i,j,ispec,i_sls) = Unp1
+ endif
if(stage_time_scheme == 6) 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_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(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_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
+ 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- &
+ 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
-
+ 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
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-01-31 04:55:42 UTC (rev 19529)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-01-31 12:44:59 UTC (rev 19530)
@@ -969,6 +969,8 @@
POWER_DISPLAY_COLOR,PERFORM_CUTHILL_MCKEE,SU_FORMAT,USER_T0, time_stepping_scheme, &
ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
if(nproc_read_from_database < 1) stop 'should have nproc_read_from_database >= 1'
+ if(SIMULATION_TYPE == 2 .and.(time_stepping_scheme == 2 .or. time_stepping_scheme == 3)) &
+ stop 'RK and LDDRK time scheme not supported for adjoint inversion'
if(nproc /= nproc_read_from_database) stop 'must always have nproc == nproc_read_from_database'
!! DK DK Dec 2011: add a crack manually
npgeo_ori = npgeo
@@ -3755,11 +3757,19 @@
! solid/porous edge detection
! the two elements forming an edge are already known (computed in meshfem2D),
! the common nodes forming the edge are computed here
- if(coupled_elastic_poro) then
+if(ATTENUATION_PORO_FLUID_PART .and. any_poroelastic .and. &
+(time_stepping_scheme == 2.or. time_stepping_scheme == 3)) &
+ stop 'RK and LDDRK time scheme not supported poroelastic simulations with attenuation'
+
+if(coupled_elastic_poro) then
+
if(ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART) &
stop 'Attenuation not supported for mixed elastic/poroelastic simulations'
+ if(time_stepping_scheme == 2.or. time_stepping_scheme == 3) &
+ stop 'RK and LDDRK time scheme not supported for mixed elastic/poroelastic simulations'
+
if ( myrank == 0 ) then
print *
print *,'Mixed elastic/poroelastic simulation'
@@ -6607,223 +6617,223 @@
endif
- if(time_stepping_scheme == 2)then
+! 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)
+! 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)
+! 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)
+! 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
+! 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
+! 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
+! 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
+! 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 .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)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)
+! 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
+! 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(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))
+! 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(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))
+! 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
+! 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(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
+! 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 .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)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)
+! 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
+! 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(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))
+! 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(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))
+! 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
- 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(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
+! 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)
+! 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
! 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
+ ! 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 .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)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)
+ ! 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
+ ! 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(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))
- 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(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))
- 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
+ ! 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)
+ ! 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 .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)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)
+ ! 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
+ ! 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(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))
+ ! 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))
- 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))
+ ! endif
- 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
+ ! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(SIMULATION_TYPE == 2) then
More information about the CIG-COMMITS
mailing list