[cig-commits] r20784 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Thu Sep 27 02:02:02 PDT 2012
Author: xie.zhinan
Date: 2012-09-27 02:02:01 -0700 (Thu, 27 Sep 2012)
New Revision: 20784
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
add mpi support for elastic simulation with LDDRK and add ADE-CPML support with LDDRK
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-09-26 23:15:02 UTC (rev 20783)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-09-27 09:02:01 UTC (rev 20784)
@@ -60,12 +60,14 @@
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,&
- e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
+ e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,nadj_rec_local, &
is_PML,nspec_PML,spec_to_PML,which_PML_elem, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
rmemory_displ_elastic,rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
+ rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
+ rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
PML_BOUNDARY_CONDITIONS)
! compute forces for the elastic elements
@@ -221,9 +223,12 @@
logical :: PML_BOUNDARY_CONDITIONS
real(kind=CUSTOM_REAL), dimension(2,3,NGLLX,NGLLZ,nspec_PML) :: rmemory_displ_elastic
+ real(kind=CUSTOM_REAL), dimension(2,3,NGLLX,NGLLZ,nspec_PML) :: rmemory_displ_elastic_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
+ rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLZ,nspec_PML) :: accel_elastic_PML
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) ::PML_dux_dxl,PML_dux_dzl,PML_duz_dxl,PML_duz_dzl,&
@@ -417,6 +422,9 @@
!---------------------- A8--------------------------
A8 = - d_x_store(i,j,ispec_PML) / (k_x_store(i,j,ispec_PML) ** 2)
bb = d_x_store(i,j,ispec_PML) / k_x_store(i,j,ispec_PML) + alpha_x_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
+
coef0 = exp(-bb * deltat)
if ( abs(bb) > 1.d-3 ) then
coef1 = (1.d0 - exp(-bb * deltat / 2.d0)) / bb
@@ -433,12 +441,33 @@
rmemory_duz_dx(i,j,ispec_PML) = coef0 * rmemory_duz_dx(i,j,ispec_PML) &
+ PML_duz_dxl_new(i,j,ispec_PML) * coef1 + PML_duz_dxl(i,j,ispec_PML) * coef2
+ endif
+
+ if(stage_time_scheme == 6) then
+
+ rmemory_dux_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_dux_dx(i,j,ispec_PML) + PML_dux_dxl(i,j,ispec_PML))
+ rmemory_dux_dx(i,j,ispec_PML) = rmemory_dux_dx(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML)
+
+ rmemory_duz_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_duz_dx(i,j,ispec_PML) + PML_duz_dxl(i,j,ispec_PML))
+ rmemory_duz_dx(i,j,ispec_PML) = rmemory_duz_dx(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML)
+
+ end if
+
dux_dxl = PML_dux_dxl(i,j,ispec_PML) + A8 * rmemory_dux_dx(i,j,ispec_PML)
duz_dxl = PML_duz_dxl(i,j,ispec_PML) + A8 * rmemory_duz_dx(i,j,ispec_PML)
+
+
!---------------------- A5--------------------------
A5 = d_x_store(i,j,ispec_PML)
bb = alpha_x_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
+
coef0 = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -457,9 +486,26 @@
rmemory_duz_dz(i,j,ispec_PML) = coef0 * rmemory_duz_dz(i,j,ispec_PML) &
+ PML_duz_dzl_new(i,j,ispec_PML) *coef1 + PML_duz_dzl(i,j,ispec_PML) * coef2
+ end if
+
+ if(stage_time_scheme == 6) then
+ rmemory_dux_dz_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dz_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_dux_dz(i,j,ispec_PML) + PML_dux_dzl(i,j,ispec_PML))
+ rmemory_dux_dz(i,j,ispec_PML) = rmemory_dux_dz(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_dux_dz_LDDRK(i,j,ispec_PML)
+
+ rmemory_duz_dz_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_duz_dz(i,j,ispec_PML) + PML_duz_dzl(i,j,ispec_PML))
+ rmemory_duz_dz(i,j,ispec_PML) = rmemory_duz_dz(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML)
+ end if
+
dux_dzl = PML_dux_dzl(i,j,ispec_PML) + A5 * rmemory_dux_dz(i,j,ispec_PML)
duz_dzl = PML_duz_dzl(i,j,ispec_PML) + A5 * rmemory_duz_dz(i,j,ispec_PML)
+
+
+
!------------------------------------------------------------------------------
!---------------------------- CORNER ------------------------------------------
!------------------------------------------------------------------------------
@@ -472,6 +518,9 @@
/ (k_x_store(i,j,ispec_PML)**2)
bb = d_x_store(i,j,ispec_PML) / k_x_store(i,j,ispec_PML) + alpha_x_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
+
coef0 = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -490,15 +539,34 @@
rmemory_duz_dx(i,j,ispec_PML) = coef0*rmemory_duz_dx(i,j,ispec_PML) &
+ PML_duz_dxl_new(i,j,ispec_PML) * coef1 + PML_duz_dxl(i,j,ispec_PML) * coef2
+ end if
+
+ if(stage_time_scheme == 6) then
+ rmemory_dux_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_dux_dx(i,j,ispec_PML) + PML_dux_dxl(i,j,ispec_PML))
+ rmemory_dux_dx(i,j,ispec_PML) = rmemory_dux_dx(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML)
+
+ rmemory_duz_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_duz_dx(i,j,ispec_PML) + PML_duz_dxl(i,j,ispec_PML))
+ rmemory_duz_dx(i,j,ispec_PML) = rmemory_duz_dx(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML)
+ end if
+
dux_dxl = PML_dux_dxl(i,j,ispec_PML) + A8 * rmemory_dux_dx(i,j,ispec_PML)
duz_dxl = PML_duz_dxl(i,j,ispec_PML) + A8 * rmemory_duz_dx(i,j,ispec_PML)
+
+
!---------------------------- A5 ----------------------------
A5 =(k_z_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML) &
- k_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML)) &
/ (k_z_store(i,j,ispec_PML)**2)
bb = d_z_store(i,j,ispec_PML) / k_z_store(i,j,ispec_PML) + alpha_z_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
+
coef0 = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -516,11 +584,25 @@
!! DK DK new from Wang eq (21)
rmemory_duz_dz(i,j,ispec_PML) = coef0 * rmemory_duz_dz(i,j,ispec_PML) &
+ PML_duz_dzl_new(i,j,ispec_PML) *coef1 + PML_duz_dzl(i,j,ispec_PML) * coef2
+ end if
+ if(stage_time_scheme == 6) then
+ rmemory_dux_dz_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dz_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_dux_dz(i,j,ispec_PML) + PML_dux_dzl(i,j,ispec_PML))
+ rmemory_dux_dz(i,j,ispec_PML) = rmemory_dux_dz(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_dux_dz_LDDRK(i,j,ispec_PML)
+ rmemory_duz_dz_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_duz_dz(i,j,ispec_PML) + PML_duz_dzl(i,j,ispec_PML))
+ rmemory_duz_dz(i,j,ispec_PML) = rmemory_duz_dz(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML)
+ end if
+
dux_dzl = PML_dux_dzl(i,j,ispec_PML) + A5 * rmemory_dux_dz(i,j,ispec_PML)
duz_dzl = PML_duz_dzl(i,j,ispec_PML) + A5 * rmemory_duz_dz(i,j,ispec_PML)
+
+
else
!------------------------------------------------------------------------------
!---------------------------- TOP & BOTTOM ------------------------------------
@@ -528,6 +610,8 @@
!---------------------- A7 --------------------------
A7 = d_z_store(i,j,ispec_PML)
bb = alpha_z_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -538,6 +622,7 @@
coef2 = deltat / 2.0d0
end if
+
!! DK DK new from Wang eq (21)
rmemory_dux_dx(i,j,ispec_PML) = coef0*rmemory_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -545,13 +630,30 @@
!! DK DK new from Wang eq (21)
rmemory_duz_dx(i,j,ispec_PML) = coef0*rmemory_duz_dx(i,j,ispec_PML) &
+ PML_duz_dxl_new(i,j,ispec_PML) * coef1 + PML_duz_dxl(i,j,ispec_PML) * coef2
+ end if
+ if(stage_time_scheme == 6) then
+ rmemory_dux_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_dux_dx(i,j,ispec_PML) + PML_dux_dxl(i,j,ispec_PML))
+ rmemory_dux_dx(i,j,ispec_PML) = rmemory_dux_dx(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML)
+
+ rmemory_duz_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_duz_dx(i,j,ispec_PML) + PML_duz_dxl(i,j,ispec_PML))
+ rmemory_duz_dx(i,j,ispec_PML) = rmemory_duz_dx(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML)
+ end if
+
dux_dxl = dux_dxl + A7 * rmemory_dux_dx(i,j,ispec_PML)
duz_dxl = duz_dxl + A7 * rmemory_duz_dx(i,j,ispec_PML)
+
+
!---------------------- A6 --------------------------
A6 = - d_z_store(i,j,ispec_PML) / ( k_z_store(i,j,ispec_PML) ** 2 )
bb = d_z_store(i,j,ispec_PML) / k_z_store(i,j,ispec_PML) + alpha_z_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
coef0 = exp(-bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -569,10 +671,25 @@
!! DK DK new from Wang eq (21)
rmemory_duz_dz(i,j,ispec_PML) = coef0 * rmemory_duz_dz(i,j,ispec_PML) &
+ PML_duz_dzl_new(i,j,ispec_PML) *coef1 + PML_duz_dzl(i,j,ispec_PML) * coef2
+ end if
+ if(stage_time_scheme == 6) then
+ rmemory_dux_dz_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dz_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_dux_dz(i,j,ispec_PML) + PML_dux_dzl(i,j,ispec_PML))
+ rmemory_dux_dz(i,j,ispec_PML) = rmemory_dux_dz(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_dux_dz_LDDRK(i,j,ispec_PML)
+
+ rmemory_duz_dz_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_duz_dz(i,j,ispec_PML) + PML_duz_dzl(i,j,ispec_PML))
+ rmemory_duz_dz(i,j,ispec_PML) = rmemory_duz_dz(i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML)
+ end if
+
dux_dzl = dux_dzl + A6 * rmemory_dux_dz(i,j,ispec_PML)
duz_dzl = duz_dzl + A6 * rmemory_duz_dz(i,j,ispec_PML)
+
+
endif
endif ! PML_BOUNDARY_CONDITIONS
@@ -767,6 +884,7 @@
.not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
bb = alpha_x_store(i,j,ispec_PML)
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -784,6 +902,7 @@
rmemory_displ_elastic(1,3,i,j,ispec_PML)=coef0 * rmemory_displ_elastic(1,3,i,j,ispec_PML) &
+ displ_elastic_new(3,iglob) * coef1 + displ_elastic(3,iglob) * coef2
rmemory_displ_elastic(2,3,i,j,ispec_PML)=0.0
+ end if
!------------------------------------------------------------------------------
!-------------------------------- CORNER --------------------------------------
@@ -793,6 +912,8 @@
!------------------------------------------------------------
bb = alpha_x_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -809,9 +930,12 @@
rmemory_displ_elastic(1,3,i,j,ispec_PML)= &
coef0 * rmemory_displ_elastic(1,3,i,j,ispec_PML) &
+ displ_elastic_new(3,iglob) * coef1 + displ_elastic(3,iglob) * coef2
+ end if
!------------------------------------------------------------
bb = alpha_z_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs(bb) > 1.d-3 ) then
@@ -828,6 +952,7 @@
rmemory_displ_elastic(2,3,i,j,ispec_PML)=coef0 * rmemory_displ_elastic(2,3,i,j,ispec_PML) &
+ displ_elastic_new(3,iglob)*(it+0.5)*deltat * coef1 &
+ displ_elastic(3,iglob)*(it-0.5)*deltat * coef2
+ end if
else
!------------------------------------------------------------------------------
@@ -836,6 +961,8 @@
!------------------------------------------------------------
bb = alpha_z_store(i,j,ispec_PML)
+
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs( bb ) > 1.d-3) then
@@ -853,6 +980,7 @@
rmemory_displ_elastic(1,3,i,j,ispec_PML)=0.d0
rmemory_displ_elastic(2,3,i,j,ispec_PML)=coef0 * rmemory_displ_elastic(2,3,i,j,ispec_PML) &
+ displ_elastic_new(3,iglob) * coef1 + displ_elastic(3,iglob) * coef2
+ end if
endif
@@ -866,6 +994,27 @@
A3 = d_x_store(i,j,ispec_PML) * alpha_x_store(i,j,ispec_PML) ** 2
A4 = 0.d0
+ if(stage_time_scheme == 6) then
+
+ bb = alpha_x_store(i,j,ispec_PML)
+
+ rmemory_displ_elastic_LDDRK(1,1,i,j,ispec_PML) = &
+ alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,1,i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_displ_elastic(1,1,i,j,ispec_PML) + displ_elastic(1,iglob))
+ rmemory_displ_elastic(1,1,i,j,ispec_PML) = rmemory_displ_elastic(1,1,i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,1,i,j,ispec_PML)
+ rmemory_displ_elastic(2,1,i,j,ispec_PML) =0.d0
+
+ rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML) = &
+ alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_displ_elastic(1,3,i,j,ispec_PML) + displ_elastic(3,iglob))
+
+ rmemory_displ_elastic(1,3,i,j,ispec_PML) = rmemory_displ_elastic(1,3,i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML)
+ rmemory_displ_elastic(2,3,i,j,ispec_PML) =0.d0
+
+ end if
+
accel_elastic_PML(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( &
A0 * displ_elastic(1,iglob) + &
@@ -876,7 +1025,7 @@
accel_elastic_PML(3,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( &
A0 * displ_elastic(3,iglob) + &
- A1 *veloc_elastic(3,iglob) + &
+ A1 * veloc_elastic(3,iglob) + &
A3 * rmemory_displ_elastic(1,3,i,j,ispec_PML) + &
A4 * rmemory_displ_elastic(2,3,i,j,ispec_PML) &
)
@@ -899,6 +1048,43 @@
(it+0.5)*deltat*alpha_x_store(i,j,ispec_PML)**2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
A4 = -alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
+ if(stage_time_scheme == 6) then
+ A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
+ d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
+ -2.d0 * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
+ A4 = alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
+ end if
+
+ if(stage_time_scheme == 6) then
+
+ bb = alpha_z_store(i,j,ispec_PML)
+
+ rmemory_displ_elastic_LDDRK(1,1,i,j,ispec_PML) = &
+ alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,1,i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_displ_elastic(1,1,i,j,ispec_PML) + displ_elastic(1,iglob))
+ rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML) = &
+ alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_displ_elastic(2,1,i,j,ispec_PML) + rmemory_displ_elastic(1,1,i,j,ispec_PML))
+
+ rmemory_displ_elastic(1,1,i,j,ispec_PML) = rmemory_displ_elastic(1,1,i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,1,i,j,ispec_PML)
+ rmemory_displ_elastic(2,1,i,j,ispec_PML) = rmemory_displ_elastic(2,1,i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML)
+
+ rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML) = &
+ alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_displ_elastic(1,3,i,j,ispec_PML) + displ_elastic(3,iglob))
+ rmemory_displ_elastic_LDDRK(2,3,i,j,ispec_PML) = &
+ alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,3,i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_displ_elastic(2,3,i,j,ispec_PML) + rmemory_displ_elastic(1,3,i,j,ispec_PML))
+
+ rmemory_displ_elastic(1,3,i,j,ispec_PML) = rmemory_displ_elastic(1,3,i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML)
+ rmemory_displ_elastic(2,3,i,j,ispec_PML) = rmemory_displ_elastic(2,3,i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,3,i,j,ispec_PML)
+
+ end if
+
accel_elastic_PML(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( &
A0 * displ_elastic(1,iglob) + &
@@ -914,6 +1100,8 @@
A4 * rmemory_displ_elastic(2,3,i,j,ispec_PML) &
)
+
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
else
@@ -923,6 +1111,26 @@
A3 = 0.d0
A4 = d_z_store(i,j,ispec_PML) * alpha_z_store(i,j,ispec_PML) ** 2
+ if(stage_time_scheme == 6) then
+
+ bb = alpha_z_store(i,j,ispec_PML)
+
+ rmemory_displ_elastic(1,1,i,j,ispec_PML) =0.d0
+ rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML) = &
+ alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_displ_elastic(2,1,i,j,ispec_PML) + displ_elastic(1,iglob))
+ rmemory_displ_elastic(2,1,i,j,ispec_PML) = rmemory_displ_elastic(2,1,i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML)
+
+ rmemory_displ_elastic(1,3,i,j,ispec_PML) =0.d0
+ rmemory_displ_elastic_LDDRK(2,3,i,j,ispec_PML) = &
+ alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,3,i,j,ispec_PML) &
+ + deltat * (-bb * rmemory_displ_elastic(2,3,i,j,ispec_PML) + displ_elastic(3,iglob))
+ rmemory_displ_elastic(2,3,i,j,ispec_PML) = rmemory_displ_elastic(2,3,i,j,ispec_PML) + &
+ beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,3,i,j,ispec_PML)
+
+ end if
+
accel_elastic_PML(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( &
A0 * displ_elastic(1,iglob) + &
@@ -937,6 +1145,9 @@
A3 * rmemory_displ_elastic(1,3,i,j,ispec_PML) + &
A4 * rmemory_displ_elastic(2,3,i,j,ispec_PML) &
)
+
+
+
endif
enddo
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-09-26 23:15:02 UTC (rev 20783)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-09-27 09:02:01 UTC (rev 20784)
@@ -467,7 +467,8 @@
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
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: veloc_elastic_LDDRK,displ_elastic_LDDRK,&
+ veloc_elastic_LDDRK_temp
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: accel_elastic_rk,veloc_elastic_rk
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: veloc_elastic_initial_rk,displ_elastic_initial_rk
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_elastic_adj_coupling
@@ -1005,10 +1006,14 @@
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
rmemory_dux_dx,rmemory_duz_dx,rmemory_dux_dz,rmemory_duz_dz
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
+ rmemory_dux_dx_LDDRK,rmemory_duz_dx_LDDRK,rmemory_dux_dz_LDDRK,rmemory_duz_dz_LDDRK
+
integer, dimension(:), allocatable :: spec_to_PML,icorner_iglob
logical, dimension(:,:), allocatable :: which_PML_elem
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic_LDDRK
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rmemory_potential_acoustic
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
@@ -2595,6 +2600,7 @@
if(time_stepping_scheme==2)then
allocate(displ_elastic_LDDRK(3,nglob_elastic))
allocate(veloc_elastic_LDDRK(3,nglob_elastic))
+ allocate(veloc_elastic_LDDRK_temp(3,nglob_elastic))
endif
if(time_stepping_scheme == 3)then
@@ -2909,12 +2915,34 @@
allocate(rmemory_duz_dz(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_duz_dz'
+ if(time_stepping_scheme == 2)then
+ allocate(rmemory_displ_elastic_LDDRK(2,3,NGLLX,NGLLZ,nspec_PML),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_displ_elastic'
+ allocate(rmemory_dux_dx_LDDRK(NGLLX,NGLLZ,nspec_PML),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_dux_dx'
+ allocate(rmemory_dux_dz_LDDRK(NGLLX,NGLLZ,nspec_PML),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_dux_dz'
+ allocate(rmemory_duz_dx_LDDRK(NGLLX,NGLLZ,nspec_PML),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_duz_dx'
+ allocate(rmemory_duz_dz_LDDRK(NGLLX,NGLLZ,nspec_PML),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_duz_dz'
+ endif
+
+
rmemory_displ_elastic(:,:,:,:,:) = ZERO
rmemory_dux_dx(:,:,:) = ZERO
rmemory_dux_dz(:,:,:) = ZERO
rmemory_duz_dx(:,:,:) = ZERO
rmemory_duz_dz(:,:,:) = ZERO
+ if(time_stepping_scheme == 2)then
+ rmemory_displ_elastic_LDDRK(:,:,:,:,:) = ZERO
+ rmemory_dux_dx_LDDRK(:,:,:) = ZERO
+ rmemory_dux_dz_LDDRK(:,:,:) = ZERO
+ rmemory_duz_dx_LDDRK(:,:,:) = ZERO
+ rmemory_duz_dz_LDDRK(:,:,:) = ZERO
+ endif
+
else
allocate(rmemory_displ_elastic(1,1,1,1,1))
@@ -2923,6 +2951,14 @@
allocate(rmemory_duz_dx(1,1,1))
allocate(rmemory_duz_dz(1,1,1))
+ if(time_stepping_scheme == 2)then
+ allocate(rmemory_displ_elastic_LDDRK(1,1,1,1,1))
+ allocate(rmemory_dux_dx_LDDRK(1,1,1))
+ allocate(rmemory_dux_dz_LDDRK(1,1,1))
+ allocate(rmemory_duz_dx_LDDRK(1,1,1))
+ allocate(rmemory_duz_dz_LDDRK(1,1,1))
+ endif
+
end if
if (any_acoustic .and. nspec_PML>0) then
@@ -3361,6 +3397,7 @@
if(time_stepping_scheme == 2) then
displ_elastic_LDDRK = 0._CUSTOM_REAL
veloc_elastic_LDDRK = 0._CUSTOM_REAL
+ veloc_elastic_LDDRK_temp = 0._CUSTOM_REAL
endif
if(time_stepping_scheme == 3)then
@@ -5405,12 +5442,14 @@
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, &
- e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
+ e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,max(1,nadj_rec_local), &
is_PML,nspec_PML,spec_to_PML,which_PML_elem, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
- rmemory_displ_elastic, rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz,&
+ rmemory_displ_elastic,rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
+ rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
+ rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
PML_BOUNDARY_CONDITIONS)
if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
@@ -5970,6 +6009,23 @@
! assembling accel_elastic for elastic elements
#ifdef USE_MPI
+
+ if(time_stepping_scheme == 2)then
+ if(i_stage==1 .and. it == 1)then
+ veloc_elastic_LDDRK_temp = veloc_elastic
+ if (nproc > 1 .and. any_elastic .and. ninterface_elastic > 0) then
+ call assemble_MPI_vector_el(veloc_elastic,nglob, &
+ ninterface, ninterface_elastic,inum_interfaces_elastic, &
+ max_interface_size, max_ibool_interfaces_size_el,&
+ ibool_interfaces_elastic, nibool_interfaces_elastic, &
+ tab_requests_send_recv_elastic,buffer_send_faces_vector_el, &
+ buffer_recv_faces_vector_el, my_neighbours)
+ endif
+ endif
+ endif
+
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+
if (nproc > 1 .and. any_elastic .and. ninterface_elastic > 0) then
call assemble_MPI_vector_el(accel_elastic,nglob, &
ninterface, ninterface_elastic,inum_interfaces_elastic, &
@@ -5979,6 +6035,7 @@
buffer_recv_faces_vector_el, my_neighbours)
endif
+
if (nproc > 1 .and. any_elastic .and. ninterface_elastic > 0 .and. SIMULATION_TYPE == 2) then
call assemble_MPI_vector_el(b_accel_elastic,nglob, &
ninterface, ninterface_elastic,inum_interfaces_elastic, &
@@ -6006,16 +6063,21 @@
endif
if(time_stepping_scheme == 2)then
+
accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic_one
accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic_one
accel_elastic(3,:) = accel_elastic(3,:) * rmass_inverse_elastic_three
- veloc_elastic_LDDRK = alpha_LDDRK(i_stage) * veloc_elastic_LDDRK + deltat * accel_elastic
- displ_elastic_LDDRK = alpha_LDDRK(i_stage) * displ_elastic_LDDRK + deltat * veloc_elastic
+ veloc_elastic_LDDRK = alpha_LDDRK(i_stage) * veloc_elastic_LDDRK + deltat * accel_elastic
+ displ_elastic_LDDRK = alpha_LDDRK(i_stage) * displ_elastic_LDDRK + deltat * veloc_elastic
+ if(i_stage==1 .and. it == 1)then
+ veloc_elastic_LDDRK_temp = veloc_elastic_LDDRK_temp + beta_LDDRK(i_stage) * veloc_elastic_LDDRK
+ veloc_elastic = veloc_elastic_LDDRK_temp
+ else
+ veloc_elastic = veloc_elastic + beta_LDDRK(i_stage) * veloc_elastic_LDDRK
+ endif
+ displ_elastic = displ_elastic + beta_LDDRK(i_stage) * displ_elastic_LDDRK
- veloc_elastic = veloc_elastic + beta_LDDRK(i_stage) * veloc_elastic_LDDRK
- displ_elastic = displ_elastic + beta_LDDRK(i_stage) * displ_elastic_LDDRK
-
endif
if(time_stepping_scheme == 3)then
More information about the CIG-COMMITS
mailing list