[cig-commits] r21364 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Wed Feb 13 11:04:17 PST 2013
Author: xie.zhinan
Date: 2013-02-13 11:04:16 -0800 (Wed, 13 Feb 2013)
New Revision: 21364
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
add pml support for elastic kernel computation
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-02-13 17:41:56 UTC (rev 21363)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-02-13 19:04:16 UTC (rev 21364)
@@ -1430,10 +1430,10 @@
enddo ! end of loop over all spectral elements
- !
+ !
!--- absorbing boundaries
!
- if(anyabs .and. .not. PML_BOUNDARY_CONDITIONS) then
+ if(anyabs .and. .not. PML_BOUNDARY_CONDITIONS .and. backward_simulation) then
count_left=1
count_right=1
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-02-13 17:41:56 UTC (rev 21363)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-02-13 19:04:16 UTC (rev 21364)
@@ -44,7 +44,8 @@
subroutine pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
nspec_PML,is_PML,which_PML_elem,spec_to_PML, &
icorner_iglob,NELEM_PML_THICKNESS,&
- read_external_mesh,region_CPML)
+ read_external_mesh,region_CPML,&
+ SIMULATION_TYPE,PML_interior_interface,nglob_interface,SAVE_FORWARD)
implicit none
@@ -54,6 +55,9 @@
include 'mpif.h'
#endif
+ integer :: SIMULATION_TYPE,nglob_interface
+ logical, dimension(4,nspec) :: PML_interior_interface
+
integer :: nspec,nglob,nelemabs,nspec_PML,NELEM_PML_THICKNESS
logical :: anyabs
@@ -71,6 +75,7 @@
!! DK DK for CPML_element_file
logical :: read_external_mesh
integer, dimension(nspec) :: region_CPML
+ logical :: SAVE_FORWARD
!!!detection of PML elements
@@ -152,10 +157,75 @@
end do !end nelem_thickness loop
+ if(SIMULATION_TYPE == 2 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))then
+
+ do i_coef=NELEM_PML_THICKNESS,NELEM_PML_THICKNESS+1
+ do ispec=1,nspec
+ if (.not. which_PML_elem(ibound,ispec)) then
+ do j=1,NGLLZ,NGLLZ-1
+ do i=1,NGLLX,NGLLX-1
+ iglob=ibool(i,j,ispec)
+ do k=1,ncorner
+ if (iglob==icorner_iglob(k)) then
+ PML_interior_interface(ibound,ispec) = .true.
+ end if
+ end do
+ end do
+ end do
+ end if
+ end do
+
+ end do !end nelem_thickness loop
+
+ endif !end of SIMULATION_TYPE == 2
+
write(IOUT,*) "number of PML spectral elements on side ", ibound,":", nspec_PML
enddo ! end loop on the 4 boundaries
+ if(SIMULATION_TYPE == 2 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))then
+ nglob_interface = 0
+ do ispec = 1,nspec
+ if(PML_interior_interface(IBOTTOM,ispec) &
+ .and. (.not. PML_interior_interface(IRIGHT,ispec)) &
+ .and. (.not. PML_interior_interface(ILEFT,ispec)) &
+ .and. (.not. which_PML_elem(IRIGHT,ispec)) &
+ .and. (.not. which_PML_elem(ILEFT,ispec)))then
+ nglob_interface = nglob_interface + 5
+ elseif(PML_interior_interface(ITOP,ispec) &
+ .and. (.not. PML_interior_interface(IRIGHT,ispec)) &
+ .and. (.not. PML_interior_interface(ILEFT,ispec)) &
+ .and. (.not. which_PML_elem(IRIGHT,ispec)) &
+ .and. (.not. which_PML_elem(ILEFT,ispec)))then
+ nglob_interface = nglob_interface + 5
+ elseif(PML_interior_interface(IRIGHT,ispec) &
+ .and. (.not. PML_interior_interface(IBOTTOM,ispec)) &
+ .and. (.not. PML_interior_interface(ITOP,ispec)) &
+ .and. (.not. which_PML_elem(IBOTTOM,ispec)) &
+ .and. (.not. which_PML_elem(ITOP,ispec)))then
+ nglob_interface = nglob_interface + 5
+ elseif(PML_interior_interface(ILEFT,ispec) &
+ .and. (.not. PML_interior_interface(IBOTTOM,ispec)) &
+ .and. (.not. PML_interior_interface(ITOP,ispec)) &
+ .and. (.not. which_PML_elem(IBOTTOM,ispec)) &
+ .and. (.not. which_PML_elem(ITOP,ispec)))then
+ nglob_interface = nglob_interface + 5
+ elseif(PML_interior_interface(ILEFT,ispec) &
+ .and. PML_interior_interface(IBOTTOM,ispec))then
+ nglob_interface = nglob_interface + 10
+ elseif(PML_interior_interface(IRIGHT,ispec) &
+ .and. PML_interior_interface(IBOTTOM,ispec))then
+ nglob_interface = nglob_interface + 10
+ elseif(PML_interior_interface(ILEFT,ispec) &
+ .and. PML_interior_interface(ITOP,ispec))then
+ nglob_interface = nglob_interface + 10
+ elseif(PML_interior_interface(IRIGHT,ispec) &
+ .and. PML_interior_interface(ITOP,ispec))then
+ nglob_interface = nglob_interface + 10
+ endif
+ enddo
+endif
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do ispec=1,nspec
if(is_PML(ispec)) then
@@ -257,7 +327,124 @@
endif
end subroutine pml_init
+!
+!-------------------------------------------------------------------------------------------------
+!
+ subroutine determin_interface_pml_interior(nglob_interface,nspec,ibool,PML_interior_interface,&
+ which_PML_elem,point_interface)
+ implicit none
+ include 'constants.h'
+
+ integer :: nglob_interface, nspec
+ logical, dimension(4,nspec) :: PML_interior_interface
+ logical, dimension(4,nspec) :: which_PML_elem
+ integer :: ispec
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+ integer, dimension(nglob_interface) :: point_interface
+
+ nglob_interface = 0
+ do ispec = 1,nspec
+ if(PML_interior_interface(IBOTTOM,ispec) &
+ .and. (.not. PML_interior_interface(IRIGHT,ispec)) &
+ .and. (.not. PML_interior_interface(ILEFT,ispec)) &
+ .and. (.not. which_PML_elem(IRIGHT,ispec)) &
+ .and. (.not. which_PML_elem(ILEFT,ispec)))then
+ point_interface(nglob_interface + 1) = ibool(1,1,ispec)
+ point_interface(nglob_interface + 2) = ibool(2,1,ispec)
+ point_interface(nglob_interface + 3) = ibool(3,1,ispec)
+ point_interface(nglob_interface + 4) = ibool(4,1,ispec)
+ point_interface(nglob_interface + 5) = ibool(5,1,ispec)
+ nglob_interface = nglob_interface + 5
+ elseif(PML_interior_interface(ITOP,ispec) &
+ .and. (.not. PML_interior_interface(IRIGHT,ispec)) &
+ .and. (.not. PML_interior_interface(ILEFT,ispec)) &
+ .and. (.not. which_PML_elem(IRIGHT,ispec)) &
+ .and. (.not. which_PML_elem(ILEFT,ispec)))then
+ point_interface(nglob_interface + 1) = ibool(1,NGLLZ,ispec)
+ point_interface(nglob_interface + 2) = ibool(2,NGLLZ,ispec)
+ point_interface(nglob_interface + 3) = ibool(3,NGLLZ,ispec)
+ point_interface(nglob_interface + 4) = ibool(4,NGLLZ,ispec)
+ point_interface(nglob_interface + 5) = ibool(5,NGLLZ,ispec)
+ nglob_interface = nglob_interface + 5
+ elseif(PML_interior_interface(IRIGHT,ispec) &
+ .and. (.not. PML_interior_interface(IBOTTOM,ispec)) &
+ .and. (.not. PML_interior_interface(ITOP,ispec)) &
+ .and. (.not. which_PML_elem(IBOTTOM,ispec)) &
+ .and. (.not. which_PML_elem(ITOP,ispec)))then
+ point_interface(nglob_interface + 1) = ibool(NGLLX,1,ispec)
+ point_interface(nglob_interface + 2) = ibool(NGLLX,2,ispec)
+ point_interface(nglob_interface + 3) = ibool(NGLLX,3,ispec)
+ point_interface(nglob_interface + 4) = ibool(NGLLX,4,ispec)
+ point_interface(nglob_interface + 5) = ibool(NGLLX,5,ispec)
+ nglob_interface = nglob_interface + 5
+ elseif(PML_interior_interface(ILEFT,ispec) &
+ .and. (.not. PML_interior_interface(IBOTTOM,ispec)) &
+ .and. (.not. PML_interior_interface(ITOP,ispec)) &
+ .and. (.not. which_PML_elem(IBOTTOM,ispec)) &
+ .and. (.not. which_PML_elem(ITOP,ispec)))then
+ point_interface(nglob_interface + 1) = ibool(1,1,ispec)
+ point_interface(nglob_interface + 2) = ibool(1,2,ispec)
+ point_interface(nglob_interface + 3) = ibool(1,3,ispec)
+ point_interface(nglob_interface + 4) = ibool(1,4,ispec)
+ point_interface(nglob_interface + 5) = ibool(1,5,ispec)
+ nglob_interface = nglob_interface + 5
+ elseif(PML_interior_interface(ILEFT,ispec) &
+ .and. PML_interior_interface(IBOTTOM,ispec))then
+ point_interface(nglob_interface + 1) = ibool(1,1,ispec)
+ point_interface(nglob_interface + 2) = ibool(1,2,ispec)
+ point_interface(nglob_interface + 3) = ibool(1,3,ispec)
+ point_interface(nglob_interface + 4) = ibool(1,4,ispec)
+ point_interface(nglob_interface + 5) = ibool(1,5,ispec)
+ point_interface(nglob_interface + 6) = ibool(1,1,ispec)
+ point_interface(nglob_interface + 7) = ibool(2,1,ispec)
+ point_interface(nglob_interface + 8) = ibool(3,1,ispec)
+ point_interface(nglob_interface + 9) = ibool(4,1,ispec)
+ point_interface(nglob_interface + 10)= ibool(5,1,ispec)
+ nglob_interface = nglob_interface + 10
+ elseif(PML_interior_interface(IRIGHT,ispec) &
+ .and. PML_interior_interface(IBOTTOM,ispec))then
+ point_interface(nglob_interface + 1) = ibool(NGLLX,1,ispec)
+ point_interface(nglob_interface + 2) = ibool(NGLLX,2,ispec)
+ point_interface(nglob_interface + 3) = ibool(NGLLX,3,ispec)
+ point_interface(nglob_interface + 4) = ibool(NGLLX,4,ispec)
+ point_interface(nglob_interface + 5) = ibool(NGLLX,5,ispec)
+ point_interface(nglob_interface + 6) = ibool(1,1,ispec)
+ point_interface(nglob_interface + 7) = ibool(2,1,ispec)
+ point_interface(nglob_interface + 8) = ibool(3,1,ispec)
+ point_interface(nglob_interface + 9) = ibool(4,1,ispec)
+ point_interface(nglob_interface + 10)= ibool(5,1,ispec)
+ nglob_interface = nglob_interface + 10
+ elseif(PML_interior_interface(ILEFT,ispec) &
+ .and. PML_interior_interface(ITOP,ispec))then
+ point_interface(nglob_interface + 1) = ibool(1,1,ispec)
+ point_interface(nglob_interface + 2) = ibool(1,2,ispec)
+ point_interface(nglob_interface + 3) = ibool(1,3,ispec)
+ point_interface(nglob_interface + 4) = ibool(1,4,ispec)
+ point_interface(nglob_interface + 5) = ibool(1,5,ispec)
+ point_interface(nglob_interface + 6) = ibool(1,NGLLZ,ispec)
+ point_interface(nglob_interface + 7) = ibool(2,NGLLZ,ispec)
+ point_interface(nglob_interface + 8) = ibool(3,NGLLZ,ispec)
+ point_interface(nglob_interface + 9) = ibool(4,NGLLZ,ispec)
+ point_interface(nglob_interface + 10)= ibool(5,NGLLZ,ispec)
+ nglob_interface = nglob_interface + 10
+ elseif(PML_interior_interface(IRIGHT,ispec) &
+ .and. PML_interior_interface(ITOP,ispec))then
+ point_interface(nglob_interface + 1) = ibool(NGLLX,1,ispec)
+ point_interface(nglob_interface + 2) = ibool(NGLLX,2,ispec)
+ point_interface(nglob_interface + 3) = ibool(NGLLX,3,ispec)
+ point_interface(nglob_interface + 4) = ibool(NGLLX,4,ispec)
+ point_interface(nglob_interface + 5) = ibool(NGLLX,5,ispec)
+ point_interface(nglob_interface + 6) = ibool(1,NGLLZ,ispec)
+ point_interface(nglob_interface + 7) = ibool(2,NGLLZ,ispec)
+ point_interface(nglob_interface + 8) = ibool(3,NGLLZ,ispec)
+ point_interface(nglob_interface + 9) = ibool(4,NGLLZ,ispec)
+ point_interface(nglob_interface + 10)= ibool(5,NGLLZ,ispec)
+ nglob_interface = nglob_interface + 10
+ endif
+ enddo
+
+ end subroutine determin_interface_pml_interior
!
!-------------------------------------------------------------------------------------------------
!
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-02-13 17:41:56 UTC (rev 21363)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-02-13 19:04:16 UTC (rev 21364)
@@ -1028,6 +1028,14 @@
integer, dimension(:), allocatable :: spec_to_PML,icorner_iglob
logical, dimension(:,:), allocatable :: which_PML_elem
+!ZN defined for using PML in elastic simulation
+ logical, dimension(:,:), allocatable :: PML_interior_interface
+ integer, dimension(:), allocatable :: point_interface
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_displ
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_veloc
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_accel
+ integer :: nglob_interface
+
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic_LDDRK
@@ -2923,13 +2931,57 @@
allocate(which_PML_elem(4,nspec))
allocate(spec_to_PML(nspec))
+ if(SIMULATION_TYPE == 2 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))then
+ allocate(PML_interior_interface(4,nspec))
+ PML_interior_interface = .false.
+ else
+ allocate(PML_interior_interface(4,1))
+ endif
+
+
is_PML(:) = .false.
which_PML_elem(:,:) = .false.
! DK DK add support for using pml in mpi mode with external mesh
call pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
nspec_PML,is_PML,which_PML_elem,spec_to_PML, &
icorner_iglob,NELEM_PML_THICKNESS,&
- read_external_mesh,region_CPML)
+ read_external_mesh,region_CPML,&
+ SIMULATION_TYPE,PML_interior_interface,nglob_interface,SAVE_FORWARD)
+
+ if((SIMULATION_TYPE == 2 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD)) .and. PML_BOUNDARY_CONDITIONS)then
+ allocate(point_interface(nglob_interface))
+ allocate(pml_interfeace_history_displ(3,nglob_interface,NSTEP))
+ allocate(pml_interfeace_history_veloc(3,nglob_interface,NSTEP))
+ allocate(pml_interfeace_history_accel(3,nglob_interface,NSTEP))
+
+ call determin_interface_pml_interior(nglob_interface,nspec,ibool,PML_interior_interface,&
+ which_PML_elem,point_interface)
+ deallocate(PML_interior_interface)
+ write(outputname,'(a,i6.6,a)') 'pml_elastic',myrank,'.bin'
+ open(unit=71,file='OUTPUT_FILES/'//outputname,status='unknown',form='unformatted')
+ else
+ allocate(point_interface(1))
+ allocate(pml_interfeace_history_displ(3,1,1))
+ allocate(pml_interfeace_history_veloc(3,1,1))
+ allocate(pml_interfeace_history_accel(3,1,1))
+ endif
+
+ if(SIMULATION_TYPE == 2 .and. PML_BOUNDARY_CONDITIONS)then
+ do it = 1,NSTEP
+ do i = 1, nglob_interface
+ read(71)pml_interfeace_history_accel(1,i,it),&
+ pml_interfeace_history_accel(2,i,it),&
+ pml_interfeace_history_accel(3,i,it),&
+ pml_interfeace_history_veloc(1,i,it),&
+ pml_interfeace_history_veloc(2,i,it),&
+ pml_interfeace_history_veloc(3,i,it),&
+ pml_interfeace_history_displ(1,i,it),&
+ pml_interfeace_history_displ(2,i,it),&
+ pml_interfeace_history_displ(3,i,it)
+ enddo
+ enddo
+ endif
+
deallocate(which_PML_elem)
deallocate(icorner_iglob)
@@ -5596,6 +5648,17 @@
if(SIMULATION_TYPE == 2)then
+ if(PML_BOUNDARY_CONDITIONS)then
+ do i = 1, nglob_interface
+ b_veloc_elastic(1,point_interface(i)) = pml_interfeace_history_veloc(1,i,NSTEP-it+1)
+ b_veloc_elastic(2,point_interface(i)) = pml_interfeace_history_veloc(2,i,NSTEP-it+1)
+ b_veloc_elastic(3,point_interface(i)) = pml_interfeace_history_veloc(3,i,NSTEP-it+1)
+ b_displ_elastic(1,point_interface(i)) = pml_interfeace_history_displ(1,i,NSTEP-it+1)
+ b_displ_elastic(2,point_interface(i)) = pml_interfeace_history_displ(2,i,NSTEP-it+1)
+ b_displ_elastic(3,point_interface(i)) = pml_interfeace_history_displ(3,i,NSTEP-it+1)
+ enddo
+ endif
+
call compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
source_type,it,NSTEP,anyabs,assign_external_model, &
@@ -5626,8 +5689,30 @@
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
- PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.)
+!ZN PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.)
+ .false.,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.)
+ do ispec = 1,nspec
+ do i = 1, NGLLX
+ do j = 1, NGLLZ
+ if(elastic(ispec) .and. is_pml(ispec))then
+ b_accel_elastic(:,ibool(i,j,ispec)) = 0.
+ b_veloc_elastic(:,ibool(i,j,ispec)) = 0.
+ b_displ_elastic(:,ibool(i,j,ispec)) = 0.
+ endif
+ enddo
+ enddo
+ enddo
+
+ if(PML_BOUNDARY_CONDITIONS)then
+ do i = 1, nglob_interface
+ b_accel_elastic(1,point_interface(i)) = pml_interfeace_history_accel(1,i,NSTEP-it+1)
+ b_accel_elastic(2,point_interface(i)) = pml_interfeace_history_accel(2,i,NSTEP-it+1)
+ b_accel_elastic(3,point_interface(i)) = pml_interfeace_history_accel(3,i,NSTEP-it+1)
+ enddo
+ endif
+
+
call compute_forces_viscoelastic_pre_kernel(p_sv,nglob,nspec,displ_elastic,b_displ_elastic,&
mu_k,kappa_k,elastic,ibool,hprime_xx,hprime_zz,xix,xiz,gammax,gammaz,SIMULATION_TYPE)
endif
@@ -5708,6 +5793,17 @@
endif ! if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)
+ if(PML_BOUNDARY_CONDITIONS .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)then
+ do i = 1, nglob_interface
+ write(71)accel_elastic(1,point_interface(i)),accel_elastic(2,point_interface(i)),&
+ accel_elastic(3,point_interface(i)),&
+ veloc_elastic(1,point_interface(i)),veloc_elastic(2,point_interface(i)),&
+ veloc_elastic(3,point_interface(i)),&
+ displ_elastic(1,point_interface(i)),displ_elastic(2,point_interface(i)),&
+ displ_elastic(3,point_interface(i))
+ enddo
+ endif !(PML_BOUNDARY_CONDITIONS .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)
+
endif !if(any_elastic)
! *********************************************************
@@ -8752,6 +8848,8 @@
close(36)
close(37)
close(38)
+ close(71)
+
endif
if(any_poroelastic) then
close(25)
More information about the CIG-COMMITS
mailing list