[cig-commits] r22918 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Thu Oct 3 02:25:55 PDT 2013
Author: xie.zhinan
Date: 2013-10-03 02:25:55 -0700 (Thu, 03 Oct 2013)
New Revision: 22918
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
Log:
fix one error in invert_mass_matrix.F90 and clean the code of compute_forces_viscoelastic.F90
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-10-03 00:20:11 UTC (rev 22917)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-10-03 09:25:55 UTC (rev 22918)
@@ -957,7 +957,6 @@
endif ! end of test if elastic element
enddo ! end of loop over all spectral elements
-
!
!--- absorbing boundaries
!
@@ -971,6 +970,7 @@
! get elastic parameters of current spectral element
lambdal_unrelaxed_elastic = poroelastcoef(1,1,kmato(ispec))
mul_unrelaxed_elastic = poroelastcoef(2,1,kmato(ispec))
+ lambdaplus2mu_unrelaxed_elastic = lambdal_unrelaxed_elastic + 2._CUSTOM_REAL * mul_unrelaxed_elastic
rhol = density(1,kmato(ispec))
kappal = lambdal_unrelaxed_elastic + TWO*mul_unrelaxed_elastic/3._CUSTOM_REAL
cpl = sqrt((kappal + 4._CUSTOM_REAL*mul_unrelaxed_elastic/3._CUSTOM_REAL)/rhol)
@@ -991,8 +991,8 @@
call compute_Bielak_conditions(coord,iglob,nglob,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
x0_source, z0_source, A_plane, B_plane, C_plane, anglesource(1), anglesource_refl, &
c_inc, c_refl, time_offset,f0)
- traction_x_t0 = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)*dxUx + lambdal_unrelaxed_elastic*dzUz
- traction_z_t0 = mul_unrelaxed_elastic*(dxUz + dzUx)
+ traction_x_t0 = lambdaplus2mu_unrelaxed_elastic * dxUx + lambdal_unrelaxed_elastic * dzUz
+ traction_z_t0 = mul_unrelaxed_elastic * (dxUz + dzUx)
else
veloc_horiz=v0x_left(count_left)
veloc_vert=v0z_left(count_left)
@@ -1045,11 +1045,11 @@
displn = nx*displx+nz*displz
- displtx = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/(2.0*spring_position)*displn*nx + &
- mul_unrelaxed_elastic/(2.0*spring_position)*(displx-displn*nx)
+ displtx = lambdaplus2mu_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * displn * nx + &
+ mul_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * (displx-displn*nx)
displty = mul_unrelaxed_elastic*disply
- displtz = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/(2.0*spring_position)*displn*nz + &
- mul_unrelaxed_elastic/(2.0*spring_position)*(displz-displn*nz)
+ displtz = lambdaplus2mu_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * displn * nz + &
+ mul_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * (displz-displn*nz)
endif
if((SIMULATION_TYPE == 3 .and. .not. backward_simulation) .or. SIMULATION_TYPE == 1) then
@@ -1095,7 +1095,7 @@
call compute_Bielak_conditions(coord,iglob,nglob,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
x0_source, z0_source, A_plane, B_plane, C_plane, anglesource(1), anglesource_refl, &
c_inc, c_refl, time_offset,f0)
- traction_x_t0 = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)*dxUx + lambdal_unrelaxed_elastic*dzUz
+ traction_x_t0 = lambdaplus2mu_unrelaxed_elastic * dxUx + lambdal_unrelaxed_elastic * dzUz
traction_z_t0 = mul_unrelaxed_elastic*(dxUz + dzUx)
else
veloc_horiz=v0x_right(count_right)
@@ -1148,11 +1148,11 @@
displn = nx*displx+nz*displz
- displtx = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/(2.0*spring_position)*displn*nx + &
- mul_unrelaxed_elastic/(2.0*spring_position)*(displx-displn*nx)
+ displtx = lambdaplus2mu_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * displn * nx + &
+ mul_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * (displx-displn*nx)
displty = mul_unrelaxed_elastic*disply
- displtz = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/(2.0*spring_position)*displn*nz + &
- mul_unrelaxed_elastic/(2.0*spring_position)*(displz-displn*nz)
+ displtz = lambdaplus2mu_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * displn * nz + &
+ mul_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * (displz-displn*nz)
endif
if((SIMULATION_TYPE == 3 .and. .not. backward_simulation) .or. SIMULATION_TYPE == 1) then
@@ -1203,8 +1203,8 @@
call compute_Bielak_conditions(coord,iglob,nglob,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
x0_source, z0_source, A_plane, B_plane, C_plane, anglesource(1), anglesource_refl, &
c_inc, c_refl, time_offset,f0)
- traction_x_t0 = mul_unrelaxed_elastic*(dxUz + dzUx)
- traction_z_t0 = lambdal_unrelaxed_elastic*dxUx + (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)*dzUz
+ traction_x_t0 = mul_unrelaxed_elastic * (dxUz + dzUx)
+ traction_z_t0 = lambdal_unrelaxed_elastic * dxUx + lambdaplus2mu_unrelaxed_elastic * dzUz
else
veloc_horiz=v0x_bot(count_bottom)
veloc_vert=v0z_bot(count_bottom)
@@ -1262,11 +1262,11 @@
displn = nx*displx+nz*displz
- displtx = (lambdal_unrelaxed_elastic+2.0*mul_unrelaxed_elastic)/(2.0*spring_position)*displn*nx + &
- mul_unrelaxed_elastic/(2.0*spring_position)*(displx-displn*nx)
+ displtx = lambdaplus2mu_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * displn * nx + &
+ mul_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * (displx-displn*nx)
displty = mul_unrelaxed_elastic*disply
- displtz = (lambdal_unrelaxed_elastic+2.0*mul_unrelaxed_elastic)/(2.0*spring_position)*displn*nz + &
- mul_unrelaxed_elastic/(2.0*spring_position)*(displz-displn*nz)
+ displtz = lambdaplus2mu_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * displn * nz + &
+ mul_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * (displz-displn*nz)
if((codeabs(IEDGE4,ispecabs) .and. i == 1) .or. (codeabs(IEDGE2,ispecabs) .and. i == NGLLX)) then
displtx = 0._CUSTOM_REAL; displty = 0._CUSTOM_REAL; displtz = 0._CUSTOM_REAL
@@ -1321,8 +1321,8 @@
call compute_Bielak_conditions(coord,iglob,nglob,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
x0_source, z0_source, A_plane, B_plane, C_plane, anglesource(1), anglesource_refl, &
c_inc, c_refl, time_offset,f0)
- traction_x_t0 = mul_unrelaxed_elastic*(dxUz + dzUx)
- traction_z_t0 = lambdal_unrelaxed_elastic*dxUx + (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)*dzUz
+ traction_x_t0 = mul_unrelaxed_elastic * (dxUz + dzUx)
+ traction_z_t0 = lambdal_unrelaxed_elastic * dxUx + lambdaplus2mu_unrelaxed_elastic * dzUz
else
veloc_horiz = 0._CUSTOM_REAL; veloc_vert = 0._CUSTOM_REAL
traction_x_t0 = 0._CUSTOM_REAL; traction_z_t0 = 0._CUSTOM_REAL
@@ -1373,11 +1373,11 @@
displn = nx*displx+nz*displz
- displtx = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/(2.0*spring_position)*displn*nx + &
- mul_unrelaxed_elastic/(2.0*spring_position)*(displx-displn*nx)
- displty = mul_unrelaxed_elastic*disply
- displtz = (lambdal_unrelaxed_elastic+2*mul_unrelaxed_elastic)/(2.0*spring_position)*displn*nz + &
- mul_unrelaxed_elastic/(2.0*spring_position)*(displz-displn*nz)
+ displtx = lambdaplus2mu_unrelaxed_elastic / (2._CUSTOM_REAL * spring_position) * displn * nx + &
+ mul_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * (displx-displn*nx)
+ displty = mul_unrelaxed_elastic * disply
+ displtz = lambdaplus2mu_unrelaxed_elastic / (2._CUSTOM_REAL * spring_position) * displn * nz + &
+ mul_unrelaxed_elastic / (2._CUSTOM_REAL*spring_position) * (displz-displn*nz)
if((codeabs(IEDGE4,ispecabs) .and. i == 1) .or. (codeabs(IEDGE2,ispecabs) .and. i == NGLLX)) then
displtx = 0._CUSTOM_REAL; displty = 0._CUSTOM_REAL; displtz = 0._CUSTOM_REAL
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2013-10-03 00:20:11 UTC (rev 22917)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2013-10-03 09:25:55 UTC (rev 22918)
@@ -301,7 +301,7 @@
!--- DK and Zhinan Xie: one per component of the wave field i.e. one per spatial dimension.
!--- DK and Zhinan Xie: This was also suggested by Jean-Paul Ampuero in 2003.
!
- if(.not. PML_BOUNDARY_CONDITIONS .and. anyabs .and. time_stepping_scheme /= 1) then
+ if(.not. PML_BOUNDARY_CONDITIONS .and. anyabs .and. time_stepping_scheme == 1) then
count_left=1
count_right=1
count_bottom=1
More information about the CIG-COMMITS
mailing list