[cig-commits] r21881 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Wed Apr 17 09:00:20 PDT 2013
Author: xie.zhinan
Date: 2013-04-17 09:00:20 -0700 (Wed, 17 Apr 2013)
New Revision: 21881
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
Log:
fix one error in CPML implementation
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-04-17 15:07:45 UTC (rev 21880)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-04-17 16:00:20 UTC (rev 21881)
@@ -24,7 +24,7 @@
!
!=====================================================================
- subroutine create_mass_matrices(nglob,nspec,ibool,PML_CONDITIONS)
+ subroutine create_mass_matrices(nglob,nspec,ibool,PML_CONDITIONS,STACEY_ABSORBING_CONDITIONS)
! returns precomputed mass matrix in rmass array
@@ -42,6 +42,9 @@
! C-PML flag
logical :: PML_CONDITIONS
+! Stacey boundary condition flag
+ logical :: STACEY_ABSORBING_CONDITIONS
+
! local parameters
double precision :: weight
real(kind=CUSTOM_REAL) :: jacobianl
@@ -169,8 +172,8 @@
endif
! Stacey absorbing conditions (adds C*deltat/2 contribution to the mass matrices on Stacey edges)
- call create_mass_matrices_Stacey(nglob,nspec,ibool)
-
+ if(STACEY_ABSORBING_CONDITIONS)call create_mass_matrices_Stacey(nglob,nspec,ibool)
+
! ocean load mass matrix
call create_mass_matrices_ocean_load(nglob,nspec,ibool)
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-04-17 15:07:45 UTC (rev 21880)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-04-17 16:00:20 UTC (rev 21881)
@@ -40,7 +40,7 @@
my_neighbours_ext_mesh, my_nelmnts_neighbours_ext_mesh, &
my_interfaces_ext_mesh, &
ibool_interfaces_ext_mesh, nibool_interfaces_ext_mesh, &
- nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, &
+ STACEY_ABSORBING_CONDITIONS, nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, &
NSPEC2D_BOTTOM, NSPEC2D_TOP,&
ibelm_xmin, ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top, &
nodes_ibelm_xmin,nodes_ibelm_xmax,nodes_ibelm_ymin,nodes_ibelm_ymax,&
@@ -228,7 +228,7 @@
if( myrank == 0) then
write(IMAIN,*) ' ...creating mass matrix '
endif
- call create_mass_matrices(nglob_dummy,nspec,ibool,PML_CONDITIONS)
+ call create_mass_matrices(nglob_dummy,nspec,ibool,PML_CONDITIONS,STACEY_ABSORBING_CONDITIONS)
! saves the binary mesh files
call sync_all()
@@ -293,14 +293,26 @@
deallocate(xstore_dummy,ystore_dummy,zstore_dummy)
endif
- if( ELASTIC_SIMULATION ) then
- deallocate(rmass)
- deallocate(rmassx,rmassy,rmassz)
- endif
- if( ACOUSTIC_SIMULATION) then
- deallocate(rmass_acoustic)
- deallocate(rmassz_acoustic)
- endif
+ if(STACEY_ABSORBING_CONDITIONS)then
+ if( ELASTIC_SIMULATION ) then
+ deallocate(rmass)
+ deallocate(rmassx,rmassy,rmassz)
+ endif
+ if( ACOUSTIC_SIMULATION) then
+ deallocate(rmass_acoustic)
+ deallocate(rmassz_acoustic)
+ endif
+ endif
+
+ if(PML_CONDITIONS)then
+ if( ELASTIC_SIMULATION ) then
+ deallocate(rmass)
+ endif
+ if( ACOUSTIC_SIMULATION) then
+ deallocate(rmass_acoustic)
+ endif
+ endif
+
if( POROELASTIC_SIMULATION) then
deallocate(rmass_solid_poroelastic,rmass_fluid_poroelastic)
endif
@@ -315,7 +327,8 @@
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
nspec2D_bottom,nspec2D_top,ANISOTROPY)
- use generate_databases_par, only: STACEY_INSTEAD_OF_FREE_SURFACE,NGNOD,NGNOD2D
+ use generate_databases_par, only: STACEY_INSTEAD_OF_FREE_SURFACE,NGNOD,NGNOD2D,&
+ PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE
use create_regions_mesh_ext_par
implicit none
@@ -416,8 +429,9 @@
! absorbing faces
num_abs_boundary_faces = nspec2D_xmin + nspec2D_xmax + nspec2D_ymin + nspec2D_ymax + nspec2D_bottom
! adds faces of free surface if it also absorbs
- if( STACEY_INSTEAD_OF_FREE_SURFACE ) num_abs_boundary_faces = num_abs_boundary_faces + nspec2D_top
-
+ if( STACEY_INSTEAD_OF_FREE_SURFACE .or. PML_INSTEAD_OF_FREE_SURFACE)then
+ num_abs_boundary_faces = num_abs_boundary_faces + nspec2D_top
+ endif
! allocates arrays to store info for each face (assumes NGLLX=NGLLY=NGLLZ)
allocate( abs_boundary_ispec(num_abs_boundary_faces), &
abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces), &
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-04-17 15:07:45 UTC (rev 21880)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-04-17 16:00:20 UTC (rev 21881)
@@ -1642,7 +1642,7 @@
! gets damping profile
if( NPOWER >= 1 ) then
! INRIA research report section 6.1: http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
- pml_damping_profile_l = - ((NPOWER + 1) * vp * log(CPML_Rcoef) / (2.d0 * delta) * damping_factor) * dist**NPOWER
+ pml_damping_profile_l = - ((NPOWER + 1) * vp * log(CPML_Rcoef) / (2.d0 * delta)) * dist**(1.2*NPOWER)
else
call exit_mpi(myrank,'C-PML error: NPOWER must be greater than or equal to 1')
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-04-17 15:07:45 UTC (rev 21880)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-04-17 16:00:20 UTC (rev 21881)
@@ -150,21 +150,30 @@
! absorbing boundary surface
write(IOUT) num_abs_boundary_faces
- if( num_abs_boundary_faces > 0 ) then
- write(IOUT) abs_boundary_ispec
- write(IOUT) abs_boundary_ijk
- write(IOUT) abs_boundary_jacobian2Dw
- write(IOUT) abs_boundary_normal
- ! store mass matrix contributions
- if(ELASTIC_SIMULATION) then
- write(IOUT) rmassx
- write(IOUT) rmassy
- write(IOUT) rmassz
- endif
- if(ACOUSTIC_SIMULATION) then
- write(IOUT) rmassz_acoustic
- endif
- endif
+ if(PML_CONDITIONS)then
+ if( num_abs_boundary_faces > 0 ) then
+ write(IOUT) abs_boundary_ispec
+ write(IOUT) abs_boundary_ijk
+ write(IOUT) abs_boundary_jacobian2Dw
+ write(IOUT) abs_boundary_normal
+ endif
+ else
+ if( num_abs_boundary_faces > 0 ) then
+ write(IOUT) abs_boundary_ispec
+ write(IOUT) abs_boundary_ijk
+ write(IOUT) abs_boundary_jacobian2Dw
+ write(IOUT) abs_boundary_normal
+ ! store mass matrix contributions
+ if(ELASTIC_SIMULATION) then
+ write(IOUT) rmassx
+ write(IOUT) rmassy
+ write(IOUT) rmassz
+ endif
+ if(ACOUSTIC_SIMULATION) then
+ write(IOUT) rmassz_acoustic
+ endif
+ endif
+ endif
write(IOUT) nspec2D_xmin
write(IOUT) nspec2D_xmax
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-04-17 15:07:45 UTC (rev 21880)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-04-17 16:00:20 UTC (rev 21881)
@@ -69,9 +69,7 @@
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION,deltat,PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE, &
- nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ ATTENUATION,deltat,PML_CONDITIONS, &
one_minus_sum_beta,factor_common, &
one_minus_sum_beta_kappa,factor_common_kappa, &
alphaval,betaval,gammaval,&
@@ -90,7 +88,8 @@
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,ispec_is_elastic)
+ phase_ispec_inner_elastic,ispec_is_elastic,&
+ abs_boundary_ijk,num_abs_boundary_faces,phase_is_inner,ispec_is_inner,abs_boundary_ispec)
! adjoint simulations: backward/reconstructed wavefield
if( SIMULATION_TYPE == 3 ) &
@@ -101,9 +100,7 @@
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION,deltat,PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE, &
- nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ ATTENUATION,deltat,PML_CONDITIONS, &
one_minus_sum_beta,factor_common, &
one_minus_sum_beta_kappa,factor_common_kappa, &
b_alphaval,b_betaval,b_gammaval, &
@@ -122,7 +119,8 @@
b_dsdx_top,b_dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,ispec_is_elastic)
+ phase_ispec_inner_elastic,ispec_is_elastic,&
+ abs_boundary_ijk,num_abs_boundary_faces,phase_is_inner,ispec_is_inner,abs_boundary_ispec)
endif
@@ -469,13 +467,11 @@
kappastore,mustore,jacobian,ibool, &
ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
- one_minus_sum_beta_kappa,factor_common_kappa, & !ZN
+ one_minus_sum_beta_kappa,factor_common_kappa, &
alphaval,betaval,gammaval, &
- NSPEC_ATTENUATION_AB,NSPEC_ATTENUATION_AB_kappa, & !ZN
-!ZN R_xx,R_yy,R_xy,R_xz,R_yz, &
- R_trace,R_xx,R_yy,R_xy,R_xz,R_yz, & !ZN
-!ZN epsilondev_xx,epsilondev_yy,epsilondev_xy, &
- epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy, & !ZN
+ NSPEC_ATTENUATION_AB,NSPEC_ATTENUATION_AB_kappa, &
+ R_trace,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy, &
epsilondev_xz,epsilondev_yz,epsilon_trace_over_3, &
ANISOTROPY,NSPEC_ANISO, &
c11store,c12store,c13store,c14store,c15store,c16store,&
@@ -528,13 +524,11 @@
kappastore,mustore,jacobian,ibool, &
ATTENUATION,deltat, &
one_minus_sum_beta,factor_common, &
- one_minus_sum_beta_kappa,factor_common_kappa, & !ZN
+ one_minus_sum_beta_kappa,factor_common_kappa, &
b_alphaval,b_betaval,b_gammaval, &
- NSPEC_ATTENUATION_AB, NSPEC_ATTENUATION_AB_kappa, & !ZN
-!ZN b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
- b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, & !ZN
-!ZN b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
- b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, & !ZN
+ NSPEC_ATTENUATION_AB, NSPEC_ATTENUATION_AB_kappa, &
+ b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
+ b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
b_epsilondev_xz,b_epsilondev_yz,b_epsilon_trace_over_3, &
ANISOTROPY,NSPEC_ANISO, &
c11store,c12store,c13store,c14store,c15store,c16store,&
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-04-17 15:07:45 UTC (rev 21880)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-04-17 16:00:20 UTC (rev 21881)
@@ -31,9 +31,7 @@
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,&
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool, &
- ATTENUATION,deltat,PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE, &
- nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ ATTENUATION,deltat,PML_CONDITIONS, &
one_minus_sum_beta,factor_common, &
one_minus_sum_beta_kappa,factor_common_kappa, &
alphaval,betaval,gammaval, &
@@ -52,9 +50,10 @@
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,ispec_is_elastic)
+ phase_ispec_inner_elastic,ispec_is_elastic,&
+ abs_boundary_ijk,num_abs_boundary_faces,phase_is_inner,ispec_is_inner,abs_boundary_ispec)
- use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,N_SLS,SAVE_MOHO_MESH,ONE_THIRD,FOUR_THIRDS,IOUT
+ use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,N_SLS,SAVE_MOHO_MESH,ONE_THIRD,FOUR_THIRDS,IOUT,NGLLSQUARE
use pml_par
use fault_solver_dynamic, only : Kelvin_Voigt_eta
use specfem_par, only : FULL_ATTENUATION_SOLID
@@ -127,20 +126,13 @@
integer :: ispec2D_moho_top, ispec2D_moho_bot
! C-PML absorbing boundary conditions
- logical :: PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE
- integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP
- integer, dimension(nspec2D_xmin) :: ibelm_xmin
- integer, dimension(nspec2D_xmax) :: ibelm_xmax
- integer, dimension(nspec2D_ymin) :: ibelm_ymin
- integer, dimension(nspec2D_ymax) :: ibelm_ymax
- integer, dimension(NSPEC2D_BOTTOM) :: ibelm_bottom
- integer, dimension(NSPEC2D_TOP) :: ibelm_top
+ logical :: PML_CONDITIONS
- logical, dimension(NSPEC_AB) :: ispec_is_elastic
+ logical, dimension(NSPEC_AB) :: ispec_is_elastic
! local parameters
integer :: i_SLS,imodulo_N_SLS
- integer :: ispec,ispec2D,iglob,ispec_p,num_elements
+ integer :: ispec,iglob,ispec_p,num_elements
integer :: i,j,k,l
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
@@ -189,6 +181,16 @@
! local C-PML absorbing boundary conditions parameters
integer :: ispec_CPML
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! outer boundary of CPML
+ integer :: num_abs_boundary_faces
+ integer :: abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces)
+ integer :: abs_boundary_ispec(num_abs_boundary_faces)
+ integer :: igll,iface
+
imodulo_N_SLS = mod(N_SLS,3)
! choses inner/outer elements
@@ -876,172 +878,35 @@
!! DK DK or something like that
! C-PML boundary
- if( PML_CONDITIONS ) then
- ! xmin
- do ispec2D=1,nspec2D_xmin
- ispec = ibelm_xmin(ispec2D)
+ if(PML_CONDITIONS)then
+ do iface=1,num_abs_boundary_faces
- if(is_CPML(ispec) .and. ispec_is_elastic(ispec)) then
- i = 1
+ ispec = abs_boundary_ispec(iface)
- do k=1,NGLLZ
- do j=1,NGLLY
- iglob = ibool(i,j,k,ispec)
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ if( ispec_is_elastic(ispec) ) then
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
+ ! reference gll points on boundary face
+ do igll = 1,NGLLSQUARE
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
+ ! gets local indices for GLL point
+ i = abs_boundary_ijk(1,igll,iface)
+ j = abs_boundary_ijk(2,igll,iface)
+ k = abs_boundary_ijk(3,igll,iface)
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
- enddo
- enddo
- endif
- enddo
+ iglob=ibool(i,j,k,ispec)
- ! xmax
- do ispec2D=1,nspec2D_xmax
- ispec = ibelm_xmax(ispec2D)
+ accel(:,iglob) = 0.0
+ veloc(:,iglob) = 0.0
+ displ(:,iglob) = 0.0
- if(is_CPML(ispec) .and. ispec_is_elastic(ispec)) then
- i = NGLLX
+ enddo
+ endif ! ispec_is_elastic
+ endif ! ispec_is_inner
+ enddo
+ endif
- do k=1,NGLLZ
- do j=1,NGLLY
- iglob = ibool(i,j,k,ispec)
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
-
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
-
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
- enddo
- enddo
- endif
- enddo
-
- ! ymin
- do ispec2D=1,nspec2D_ymin
- ispec = ibelm_ymin(ispec2D)
-
- if(is_CPML(ispec) .and. ispec_is_elastic(ispec)) then
- j = 1
-
- do k=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
-
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
-
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
- enddo
- enddo
- endif
- enddo
-
- ! ymax
- do ispec2D=1,nspec2D_ymax
- ispec = ibelm_ymax(ispec2D)
-
- if(is_CPML(ispec) .and. ispec_is_elastic(ispec)) then
- j = NGLLY
-
- do k=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
-
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
-
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
- enddo
- enddo
- endif
- enddo
-
- ! bottom (zmin)
- do ispec2D=1,NSPEC2D_BOTTOM
- ispec = ibelm_bottom(ispec2D)
-
- if(is_CPML(ispec) .and. ispec_is_elastic(ispec)) then
- k = 1
-
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
-
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
-
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
- enddo
- enddo
- endif
- enddo
-
- ! top (zmax)
- if(PML_INSTEAD_OF_FREE_SURFACE)then
- do ispec2D=1,NSPEC2D_BOTTOM
- ispec = ibelm_top(ispec2D)
-
- if(is_CPML(ispec) .and. ispec_is_elastic(ispec)) then
- k = NGLLZ
-
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
-
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
-
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
- enddo
- enddo
- endif
- enddo
- endif
-
- endif ! if( PML_CONDITIONS )
-
end subroutine compute_forces_viscoelastic_noDev
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-04-17 15:07:45 UTC (rev 21880)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-04-17 16:00:20 UTC (rev 21881)
@@ -786,11 +786,11 @@
! A5 = temp_A5
if( ispec_is_elastic(ispec) ) then
- A3 = temp_A3 + (it+0.0) * deltat*temp_A4 !+ ((it+0.0) * deltat)**2*temp_A5
- A4 = -temp_A4 ! -2.0*(it+0.0) * deltat*temp_A5
+ A3 = temp_A3 !+ (it+0.0) * deltat*temp_A4 !+ ((it+0.0) * deltat)**2*temp_A5
+ A4 = 0.0 !-temp_A4 ! -2.0*(it+0.0) * deltat*temp_A5
else if( ispec_is_acoustic(ispec)) then
- A3 = temp_A3 + (it+0.5)*deltat*temp_A4 !+ ((it+0.5)*deltat)**2*temp_A5
- A4 = -temp_A4 !-2.0*(it+0.5)*deltat*temp_A5
+ A3 = temp_A3 !+ (it+0.5)*deltat*temp_A4 !+ ((it+0.5)*deltat)**2*temp_A5
+ A4 = 0.0 !-temp_A4 !-2.0*(it+0.5)*deltat*temp_A5
endif
A5 = 0.0 ! temp_A5
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-04-17 15:07:45 UTC (rev 21880)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-04-17 16:00:20 UTC (rev 21881)
@@ -189,9 +189,9 @@
epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY),stat=ier)
if( ier /= 0 ) stop 'error allocating array epsilondev_xx etc.'
- allocate(R_trace(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa,N_SLS),& !ZN
- epsilondev_trace(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa),stat=ier) !ZN
- if( ier /= 0 ) stop 'error allocating array R_trace etc.' !ZN
+ allocate(R_trace(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa,N_SLS),&
+ epsilondev_trace(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array R_trace etc.'
! note: needed for argument of deville routine
allocate(epsilon_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
@@ -202,9 +202,9 @@
factor_common(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array one_minus_sum_beta etc.'
- allocate(one_minus_sum_beta_kappa(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa), & !ZN
- factor_common_kappa(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa),stat=ier) !ZN
- if( ier /= 0 ) stop 'error allocating array one_minus_sum_beta_kappa etc.' !ZN
+ allocate(one_minus_sum_beta_kappa(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa), &
+ factor_common_kappa(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array one_minus_sum_beta_kappa etc.'
! reads mass matrices
read(27,iostat=ier) rmass
@@ -394,21 +394,30 @@
allocate(Tract_dsm_boundary(1,1,1,1))
endif
- if( num_abs_boundary_faces > 0 ) then
- read(27) abs_boundary_ispec
- read(27) abs_boundary_ijk
- read(27) abs_boundary_jacobian2Dw
- read(27) abs_boundary_normal
- ! reads mass matrix contributions on boundaries
- if( ELASTIC_SIMULATION ) then
- read(27) rmassx
- read(27) rmassy
- read(27) rmassz
- endif
- if( ACOUSTIC_SIMULATION ) then
- read(27) rmassz_acoustic
- endif
- endif
+ if(PML_CONDITIONS)then
+ if( num_abs_boundary_faces > 0 ) then
+ read(27) abs_boundary_ispec
+ read(27) abs_boundary_ijk
+ read(27) abs_boundary_jacobian2Dw
+ read(27) abs_boundary_normal
+ endif
+ else
+ if( num_abs_boundary_faces > 0 ) then
+ read(27) abs_boundary_ispec
+ read(27) abs_boundary_ijk
+ read(27) abs_boundary_jacobian2Dw
+ read(27) abs_boundary_normal
+ ! store mass matrix contributions
+ if(ELASTIC_SIMULATION) then
+ read(27) rmassx
+ read(27) rmassy
+ read(27) rmassz
+ endif
+ if(ACOUSTIC_SIMULATION) then
+ read(27) rmassz_acoustic
+ endif
+ endif
+ endif
read(27) nspec2D_xmin
read(27) nspec2D_xmax
@@ -791,9 +800,9 @@
allocate(b_epsilon_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
if( ier /= 0 ) stop 'error allocating array b_epsilon_trace_over_3'
- ! allocates attenuation solids for considering kappa !ZN
- allocate(b_R_trace(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa,N_SLS),& !ZN
- b_epsilondev_trace(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa),stat=ier) !ZN
+ ! allocates attenuation solids for considering kappa
+ allocate(b_R_trace(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa,N_SLS),&
+ b_epsilondev_trace(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa),stat=ier)
if( ier /= 0 ) stop 'error allocating array b_R_trace etc.'
else
More information about the CIG-COMMITS
mailing list