[cig-commits] r20372 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Jun 14 05:43:37 PDT 2012
Author: dkomati1
Date: 2012-06-14 05:43:36 -0700 (Thu, 14 Jun 2012)
New Revision: 20372
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
Zhinan Xie fixed his C-PML code
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-06-13 15:43:33 UTC (rev 20371)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-06-14 12:43:36 UTC (rev 20372)
@@ -75,7 +75,7 @@
include "constants.h"
- logical :: p_sv,PML_BOUNDARY_CONDITIONS
+ logical :: p_sv
integer :: NSOURCES, i_source
integer :: nglob,nspec,myrank,nelemabs,numat,it,NSTEP
integer, dimension(NSOURCES) :: ispec_selected_source,is_proc_source,source_type
@@ -220,6 +220,7 @@
logical, dimension(nspec) :: is_PML
integer, dimension(nspec) :: spec_to_PML
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool_PML
+ logical :: PML_BOUNDARY_CONDITIONS, anyabs_local
real(kind=CUSTOM_REAL), dimension(2,3,NGLLX,NGLLZ,nspec_PML) :: rmemory_displ_elastic,rmemory_displ_elastic_corner
real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLZ,nspec_PML) :: &
@@ -245,6 +246,7 @@
sigma_xz = 0
sigma_zy = 0
sigma_zz = 0
+ sigma_zx = 0
if( PML_BOUNDARY_CONDITIONS ) then
accel_elastic_PML = 0._CUSTOM_REAL
@@ -258,6 +260,9 @@
PML_duz_dxl_new = 0._CUSTOM_REAL
PML_duz_dzl_new = 0._CUSTOM_REAL
displ_elastic_new = displ_elastic + deltat * veloc_elastic
+ anyabs_local=.false.
+ elseif(.not.PML_BOUNDARY_CONDITIONS .and. anyabs )then
+ anyabs_local=.true.
endif
@@ -397,7 +402,6 @@
endif
endif
- iglob = ibool(i,j,ispec)
if( PML_BOUNDARY_CONDITIONS ) then
if( is_PML(ispec) ) then
@@ -1208,7 +1212,7 @@
!
!--- absorbing boundaries
!
- if(anyabs) then
+ if(anyabs_local) then
count_left=1
count_right=1
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-06-13 15:43:33 UTC (rev 20371)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-06-14 12:43:36 UTC (rev 20372)
@@ -59,7 +59,8 @@
#ifdef USE_GUENNEAU
,coord &
#endif
- )
+ ,K_x_store,K_z_store,npoin_PML,ibool_PML,is_PML,&
+ d_x_store,d_z_store,alpha_x_store,alpha_z_store,PML_BOUNDARY_CONDITIONS)
! builds the global mass matrix
@@ -88,7 +89,7 @@
! inverse mass matrices
integer :: nglob_elastic
- real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic_one,& !zhinan
+ real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic_one,&
rmass_inverse_elastic_three
integer :: nglob_acoustic
@@ -112,15 +113,21 @@
double precision, dimension(2,numat) :: density
double precision, dimension(4,3,numat) :: poroelastcoef
double precision, dimension(numat) :: porosity,tortuosity
- double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,rhoext,vsext !zhinan
+ double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,rhoext,vsext
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz !zhinan
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
! local parameters
integer :: ispec,i,j,iglob
double precision :: rhol,kappal,mul_relaxed,lambdal_relaxed
double precision :: rhol_s,rhol_f,rhol_bar,phil,tortl
+!!!!!!!!!!!!! DK DK added this
+ integer :: npoin_PML,iPML
+ real(kind=CUSTOM_REAL), dimension(npoin_PML) :: K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool_PML
+ logical, dimension(nspec) :: is_PML
+ logical :: PML_BOUNDARY_CONDITIONS, anyabs_local
!! DK DK added this for Guenneau, March 2012
#ifdef USE_GUENNEAU
double precision, dimension(NDIM,nglob_elastic), intent(in) :: coord
@@ -134,6 +141,8 @@
if(any_poroelastic) rmass_s_inverse_poroelastic(:) = 0._CUSTOM_REAL
if(any_poroelastic) rmass_w_inverse_poroelastic(:) = 0._CUSTOM_REAL
if(any_acoustic) rmass_inverse_acoustic(:) = 0._CUSTOM_REAL
+ if(PML_BOUNDARY_CONDITIONS)anyabs_local=.false.
+ if(.not.PML_BOUNDARY_CONDITIONS .and. anyabs)anyabs_local=.true.
do ispec = 1,nspec
do j = 1,NGLLZ
@@ -173,6 +182,16 @@
! for elastic medium
+
+ if(PML_BOUNDARY_CONDITIONS)then
+ if (is_PML(ispec)) then
+ iPML=ibool_PML(i,j,ispec)
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(iPML) * K_z_store(iPML)&
+ + (d_x_store(iPML)*k_z_store(iPML)+d_z_store(iPML)*k_x_store(iPML)) * deltat / 2.d0)
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ else
+
!! DK DK added this for Guenneau, March 2012
#ifdef USE_GUENNEAU
include "include_mass_Guenneau.f90"
@@ -187,8 +206,28 @@
#endif
rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ endif
else
+!! DK DK added this for Guenneau, March 2012
+#ifdef USE_GUENNEAU
+ include "include_mass_Guenneau.f90"
+#endif
+!! DK DK added this for Guenneau, March 2012
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
+
+#ifdef USE_GUENNEAU
+ endif
+#endif
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ endif
+
+
+
+
+ else
+
! for acoustic medium
rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
@@ -208,7 +247,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(anyabs) then
+ if(anyabs_local) then
count_left=1
count_right=1
count_bottom=1
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-06-13 15:43:33 UTC (rev 20371)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-06-14 12:43:36 UTC (rev 20372)
@@ -969,7 +969,7 @@
integer, dimension(:,:,:), allocatable :: ibool_PML
logical, dimension(:,:), allocatable :: which_PML_elem, which_PML_poin
- real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic_corner
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rmemory_potential_acoustic,&
@@ -979,7 +979,7 @@
logical :: anyabs_glob
integer :: nspec_PML, npoin_PML
-
+
!***********************************************************************
!
! i n i t i a l i z a t i o n p h a s e
@@ -2804,15 +2804,15 @@
allocate(K_z_store(npoin_PML))
allocate(d_x_store(npoin_PML))
allocate(d_z_store(npoin_PML))
- allocate(alpha_x_store(npoin_PML))
+ allocate(alpha_x_store(npoin_PML))
allocate(alpha_z_store(npoin_PML))
!! DK DK initialize to zero
K_x_store(:) = 0
K_z_store(:) = 0
- d_x_store(:) = 0
- d_z_store(:) = 0
- alpha_x_store(:) = 0
+ d_x_store(:) = 0
+ d_z_store(:) = 0
+ alpha_x_store(:) = 0
alpha_z_store(:) = 0
call define_PML_coefficients(nglob,nspec,is_PML,ibool,coord,&
@@ -2825,57 +2825,67 @@
!elastic PML memory variables
if (any_elastic .and. nspec_PML>0) then
- allocate(rmemory_displ_elastic(2,3,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_displ_elastic_corner(2,3,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_displ_elastic(2,3,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_displ_elastic_corner(2,3,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_dux_dx(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_dux_dz(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_duz_dx(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_duz_dz(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_dux_dx(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_dux_dz(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_duz_dx(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_duz_dz(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_dux_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_dux_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_duz_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_duz_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_dux_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_dux_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_duz_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_duz_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
+ rmemory_dux_dx(:,:,:,:) = ZERO
+ rmemory_dux_dz(:,:,:,:) = ZERO
+ rmemory_duz_dx(:,:,:,:) = ZERO
+ rmemory_duz_dz(:,:,:,:) = ZERO
+
+ rmemory_dux_dx_corner(:,:,:,:) = ZERO
+ rmemory_dux_dz_corner(:,:,:,:) = ZERO
+ rmemory_duz_dx_corner(:,:,:,:) = ZERO
+ rmemory_duz_dz_corner(:,:,:,:) = ZERO
+
else
- allocate(rmemory_displ_elastic(1,1,1,1,1))
- allocate(rmemory_displ_elastic_corner(1,1,1,1,1))
+ allocate(rmemory_displ_elastic(1,1,1,1,1))
+ allocate(rmemory_displ_elastic_corner(1,1,1,1,1))
- allocate(rmemory_dux_dx(1,1,1,1))
- allocate(rmemory_dux_dz(1,1,1,1))
- allocate(rmemory_duz_dx(1,1,1,1))
- allocate(rmemory_duz_dz(1,1,1,1))
+ allocate(rmemory_dux_dx(1,1,1,1))
+ allocate(rmemory_dux_dz(1,1,1,1))
+ allocate(rmemory_duz_dx(1,1,1,1))
+ allocate(rmemory_duz_dz(1,1,1,1))
- allocate(rmemory_dux_dx_corner(1,1,1,1))
- allocate(rmemory_dux_dz_corner(1,1,1,1))
- allocate(rmemory_duz_dx_corner(1,1,1,1))
- allocate(rmemory_duz_dz_corner(1,1,1,1))
+ allocate(rmemory_dux_dx_corner(1,1,1,1))
+ allocate(rmemory_dux_dz_corner(1,1,1,1))
+ allocate(rmemory_duz_dx_corner(1,1,1,1))
+ allocate(rmemory_duz_dz_corner(1,1,1,1))
end if
if (any_acoustic .and. nspec_PML>0) then
- allocate(rmemory_potential_acoustic(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_potential_ac_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_potential_acoustic(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_potential_ac_corner(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_acoustic_dux_dx(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_acoustic_dux_dz(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_acoustic_dux_dx(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_acoustic_dux_dz(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_acoustic_dux_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
- allocate(rmemory_acoustic_dux_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_acoustic_dux_dx_corner(2,NGLLX,NGLLZ,nspec_PML))
+ allocate(rmemory_acoustic_dux_dz_corner(2,NGLLX,NGLLZ,nspec_PML))
else
- allocate(rmemory_potential_acoustic(1,1,1,1))
- allocate(rmemory_potential_ac_corner(1,1,1,1))
+ allocate(rmemory_potential_acoustic(1,1,1,1))
+ allocate(rmemory_potential_ac_corner(1,1,1,1))
- allocate(rmemory_acoustic_dux_dx(1,1,1,1))
- allocate(rmemory_acoustic_dux_dz(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dx(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dz(1,1,1,1))
- allocate(rmemory_acoustic_dux_dx_corner(1,1,1,1))
- allocate(rmemory_acoustic_dux_dz_corner(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dx_corner(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dz_corner(1,1,1,1))
end if
@@ -2884,48 +2894,40 @@
!!!!!!!!!!!!! DK DK added this
- allocate(rmemory_dux_dx(1,1,1,1))
- allocate(rmemory_dux_dz(1,1,1,1))
- allocate(rmemory_duz_dx(1,1,1,1))
- allocate(rmemory_duz_dz(1,1,1,1))
+ allocate(rmemory_dux_dx(1,1,1,1))
+ allocate(rmemory_dux_dz(1,1,1,1))
+ allocate(rmemory_duz_dx(1,1,1,1))
+ allocate(rmemory_duz_dz(1,1,1,1))
- allocate(rmemory_dux_dx_corner(1,1,1,1))
- allocate(rmemory_dux_dz_corner(1,1,1,1))
- allocate(rmemory_duz_dx_corner(1,1,1,1))
- allocate(rmemory_duz_dz_corner(1,1,1,1))
+ allocate(rmemory_dux_dx_corner(1,1,1,1))
+ allocate(rmemory_dux_dz_corner(1,1,1,1))
+ allocate(rmemory_duz_dx_corner(1,1,1,1))
+ allocate(rmemory_duz_dz_corner(1,1,1,1))
- allocate(rmemory_displ_elastic(1,1,1,1,1))
- allocate(rmemory_displ_elastic_corner(1,1,1,1,1))
+ allocate(rmemory_displ_elastic(1,1,1,1,1))
+ allocate(rmemory_displ_elastic_corner(1,1,1,1,1))
- allocate(rmemory_acoustic_dux_dx(1,1,1,1))
- allocate(rmemory_acoustic_dux_dz(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dx(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dz(1,1,1,1))
- allocate(rmemory_acoustic_dux_dx_corner(1,1,1,1))
- allocate(rmemory_acoustic_dux_dz_corner(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dx_corner(1,1,1,1))
+ allocate(rmemory_acoustic_dux_dz_corner(1,1,1,1))
allocate(is_PML(1))
- allocate(ibool_PML(1,1,1))
- allocate(spec_to_PML(1))
- allocate(which_PML_elem(1,1))
+ allocate(ibool_PML(1,1,1))
+ allocate(spec_to_PML(1))
+ allocate(which_PML_elem(1,1))
allocate(K_x_store(1))
allocate(K_z_store(1))
allocate(d_x_store(1))
allocate(d_z_store(1))
- allocate(alpha_x_store(1))
+ allocate(alpha_x_store(1))
allocate(alpha_z_store(1))
-
+
end if ! PML_BOUNDARY_CONDITIONS
- rmemory_dux_dx(:,:,:,:) = ZERO
- rmemory_dux_dz(:,:,:,:) = ZERO
- rmemory_duz_dx(:,:,:,:) = ZERO
- rmemory_duz_dz(:,:,:,:) = ZERO
- rmemory_dux_dx_corner(:,:,:,:) = ZERO
- rmemory_dux_dz_corner(:,:,:,:) = ZERO
- rmemory_duz_dx_corner(:,:,:,:) = ZERO
- rmemory_duz_dz_corner(:,:,:,:) = ZERO
endif ! ipass == 1
@@ -2948,7 +2950,8 @@
#ifdef USE_GUENNEAU
,coord &
#endif
- )
+ ,K_x_store,K_z_store,npoin_PML,ibool_PML,is_PML,&
+ d_x_store,d_z_store,alpha_x_store,alpha_z_store,PML_BOUNDARY_CONDITIONS)
#ifdef USE_MPI
if ( nproc > 1 ) then
More information about the CIG-COMMITS
mailing list