[cig-commits] r21789 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Wed Apr 10 04:49:17 PDT 2013
Author: xie.zhinan
Date: 2013-04-10 04:49:16 -0700 (Wed, 10 Apr 2013)
New Revision: 21789
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
add pml support in adjoint simulation when using external mesh
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-04-09 23:01:09 UTC (rev 21788)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-04-10 11:49:16 UTC (rev 21789)
@@ -45,7 +45,7 @@
nspec_PML,is_PML,which_PML_elem,spec_to_PML, &
icorner_iglob,NELEM_PML_THICKNESS,&
read_external_mesh,region_CPML,&
- SIMULATION_TYPE,PML_interior_interface,nglob_interface,SAVE_FORWARD,myrank)
+ SIMULATION_TYPE,PML_interior_interface,nglob_interface,SAVE_FORWARD,myrank,mask_ibool)
implicit none
include 'constants.h'
@@ -82,6 +82,8 @@
integer :: ier
#endif
+ logical, dimension(nglob) :: mask_ibool
+
nspec_PML = 0
! detection of PML elements
@@ -180,14 +182,13 @@
end do
end if
end do
-
end do !end nelem_thickness loop
endif !end of SIMULATION_TYPE == 3
enddo ! end loop on the four boundaries
- if(SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))then
+ if(SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))then
nglob_interface = 0
do ispec = 1,nspec
if(PML_interior_interface(IBOTTOM,ispec) &
@@ -228,12 +229,10 @@
nglob_interface = nglob_interface + 10
endif
enddo
-endif
+ endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do ispec=1,nspec
+ do ispec=1,nspec
if(is_PML(ispec)) then
-
! element is in the left cpml layer
if( &
(which_PML_elem(ILEFT,ispec) .eqv. .true.) .and. &
@@ -291,15 +290,12 @@
(which_PML_elem(IBOTTOM,ispec) .eqv. .true. )) then
region_CPML(ispec) = CPML_BOTTOM_RIGHT
else
-
region_CPML(ispec) = 0
-
endif
endif
- enddo
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ enddo
- !construction of table to use less memory for absorbing coefficients
+ !construction of table to use less memory for absorbing coefficients
spec_to_PML=0
nspec_PML=0
do ispec=1,nspec
@@ -309,52 +305,85 @@
end if
enddo
- endif
+ endif !end of detection of element inside PML layer for inner mesher
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(read_external_mesh) then
- is_PML(:) = .false.
- which_PML_elem(:,:) = .false.
+ is_PML(:) = .false.
+ which_PML_elem(:,:) = .false.
+ nspec_PML = 0
+ spec_to_PML=0
+ mask_ibool(:) = .false.
+ do ispec=1,nspec
+ if(region_CPML(ispec) /= 0) then
+ nspec_PML = nspec_PML + 1
+ is_PML(ispec)=.true.
+ spec_to_PML(ispec)=nspec_PML
+ endif
+ if(SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))then
+ if(region_CPML(ispec) == 0) then
+ do i = 1, NGLLX
+ do j = 1, NGLLZ
+ iglob = ibool(i,j,ispec)
+ mask_ibool(iglob) = .true.
+ enddo
+ enddo
+ endif
+ endif
+ enddo
- nspec_PML = 0
- do ispec=1,nspec
- if(region_CPML(ispec) /= 0) then
- nspec_PML = nspec_PML + 1
- is_PML(ispec)=.true.
- spec_to_PML(ispec)=nspec_PML
- endif
- enddo
-
+ nglob_interface = 0
+ if(SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))then
+ do ispec=1,nspec
+ if(region_CPML(ispec) /= 0) then
+ do i = 1, NGLLX
+ do j = 1, NGLLZ
+ iglob = ibool(i,j,ispec)
+ if(mask_ibool(iglob))nglob_interface = nglob_interface + 1
+ enddo
+ enddo
+ endif
+ enddo
+ endif
endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef USE_MPI
call MPI_REDUCE(nspec_PML, nspec_PML_tot, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier)
#else
nspec_PML_tot = nspec_PML
#endif
+
if(myrank == 0) then
write(IOUT,*) "Total number of PML spectral elements: ", nspec_PML_tot
write(IOUT,*)
endif
- end subroutine pml_init
+end subroutine pml_init
!
!-------------------------------------------------------------------------------------------------
!
subroutine determin_interface_pml_interior(nglob_interface,nspec,ibool,PML_interior_interface,&
- which_PML_elem,point_interface)
+ which_PML_elem,point_interface,read_external_mesh,mask_ibool,region_CPML,nglob)
implicit none
include 'constants.h'
- integer :: nglob_interface, nspec
+ integer :: nglob_interface, nspec,nglob
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
+ logical :: read_external_mesh
+ logical, dimension(nglob) :: mask_ibool
+ integer, dimension(nspec) :: region_CPML
+ integer :: i,j,iglob
nglob_interface = 0
+
+ if(.not. read_external_mesh) then
do ispec = 1,nspec
if(PML_interior_interface(IBOTTOM,ispec) &
.and. (.not. PML_interior_interface(IRIGHT,ispec)) &
@@ -454,7 +483,25 @@
nglob_interface = nglob_interface + 10
endif
enddo
+ endif
+ if(read_external_mesh) then
+
+ nglob_interface = 0
+
+ do ispec=1,nspec
+ if(region_CPML(ispec) /= 0) then
+ do i = 1, NGLLX
+ do j = 1, NGLLZ
+ iglob = ibool(i,j,ispec)
+ if(mask_ibool(iglob))nglob_interface = nglob_interface + 1
+ point_interface(nglob_interface)= iglob
+ enddo
+ enddo
+ endif
+ enddo
+ endif
+
end subroutine determin_interface_pml_interior
!
!-------------------------------------------------------------------------------------------------
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-04-09 23:01:09 UTC (rev 21788)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-04-10 11:49:16 UTC (rev 21789)
@@ -2943,11 +2943,19 @@
is_PML(:) = .false.
which_PML_elem(:,:) = .false.
! DK DK add support for using pml in mpi mode with external mesh
+ if((SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD)) .and. PML_BOUNDARY_CONDITIONS)then
+ if(read_external_mesh)then
+ allocate(mask_ibool(nglob))
+ else
+ allocate(mask_ibool(1))
+ endif
+ endif
+
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,&
- SIMULATION_TYPE,PML_interior_interface,nglob_interface,SAVE_FORWARD,myrank)
+ SIMULATION_TYPE,PML_interior_interface,nglob_interface,SAVE_FORWARD,myrank,mask_ibool)
if((SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD)) .and. PML_BOUNDARY_CONDITIONS)then
allocate(point_interface(nglob_interface))
@@ -2959,8 +2967,10 @@
allocate(pml_interfeace_history_potential_dot_dot(nglob_interface,NSTEP))
call determin_interface_pml_interior(nglob_interface,nspec,ibool,PML_interior_interface,&
- which_PML_elem,point_interface)
+ which_PML_elem,point_interface,read_external_mesh,&
+ mask_ibool,region_CPML,nglob)
deallocate(PML_interior_interface)
+ deallocate(mask_ibool)
if(any_elastic .and. nglob_interface > 0)then
write(outputname,'(a,i6.6,a)') 'pml_interface_elastic',myrank,'.bin'
More information about the CIG-COMMITS
mailing list