[cig-commits] [commit] devel: added support for acoustic media for periodic boundary conditions (6c1775f)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Feb 20 13:30:40 PST 2014
Repository : ssh://geoshell/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/b662f42df413cb3e6b4870d885b2bbef66e60358...6031c1354c52525fc0645aa007ddc7dc627cc45e
>---------------------------------------------------------------
commit 6c1775ffacea3bbd5657ec588093ad53243c509b
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Thu Feb 20 22:24:54 2014 +0100
added support for acoustic media for periodic boundary conditions
>---------------------------------------------------------------
6c1775ffacea3bbd5657ec588093ad53243c509b
DATA/Par_file | 5 +--
src/specfem2D/enforce_acoustic_free_surface.f90 | 13 ++++--
src/specfem2D/get_global.f90 | 58 +++++++++----------------
src/specfem2D/specfem2D.F90 | 56 +++++++++++++++++++-----
4 files changed, 77 insertions(+), 55 deletions(-)
diff --git a/DATA/Par_file b/DATA/Par_file
index 30f1097..2188636 100644
--- a/DATA/Par_file
+++ b/DATA/Par_file
@@ -137,10 +137,7 @@ ADD_SPRING_TO_STACEY = .false.
ADD_PERIODIC_CONDITIONS = .false.
# horizontal periodicity distance for periodic conditions
-PERIODIC_horiz_dist = 0.3597d0
-
-# grid point detection tolerance for periodic conditions
-PERIODIC_DETECT_TOL = 3.3334d-6
+PERIODIC_HORIZ_DIST = 0.3597d0
#-----------------------------------------------------------------------------
# PARAMETERS FOR EXTERNAL MESHING
diff --git a/src/specfem2D/enforce_acoustic_free_surface.f90 b/src/specfem2D/enforce_acoustic_free_surface.f90
index eb6d672..cb8ee7d 100644
--- a/src/specfem2D/enforce_acoustic_free_surface.f90
+++ b/src/specfem2D/enforce_acoustic_free_surface.f90
@@ -44,7 +44,7 @@
subroutine enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,acoustic_surface, &
- ibool,nelem_acoustic_surface,nglob,nspec)
+ ibool,nelem_acoustic_surface,nglob,nspec,this_ibool_is_a_periodic_edge)
! free surface for an acoustic medium
! if acoustic, the free surface condition is a Dirichlet condition for the potential,
@@ -60,6 +60,8 @@
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+ logical, dimension(nglob) :: this_ibool_is_a_periodic_edge
+
real(kind=CUSTOM_REAL), dimension(nglob) :: &
potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
@@ -76,9 +78,12 @@
do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
iglob = ibool(i,j,ispec)
- potential_acoustic(iglob) = ZERO
- potential_dot_acoustic(iglob) = ZERO
- potential_dot_dot_acoustic(iglob) = ZERO
+ ! make sure that an acoustic free surface is not enforced on periodic edges
+ if(.not. this_ibool_is_a_periodic_edge(iglob)) then
+ potential_acoustic(iglob) = ZERO
+ potential_dot_acoustic(iglob) = ZERO
+ potential_dot_dot_acoustic(iglob) = ZERO
+ endif
enddo
enddo
diff --git a/src/specfem2D/get_global.f90 b/src/specfem2D/get_global.f90
index a898402..85ef409 100644
--- a/src/specfem2D/get_global.f90
+++ b/src/specfem2D/get_global.f90
@@ -42,29 +42,23 @@
!
!========================================================================
- subroutine get_global(nspec,nglob,ibool)
+! create a sorted version of the indirect addressing to reduce cache misses
+
+ subroutine get_global(nspec,nglob,ibool,copy_ibool_ori,integer_mask_ibool)
implicit none
include "constants.h"
integer :: nspec,nglob
- integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool,copy_ibool_ori
+ integer, dimension(nglob) :: integer_mask_ibool
! local parameters
- integer, dimension(:,:,:), allocatable :: copy_ibool_ori
- integer, dimension(:), allocatable :: mask_ibool
- integer :: inumber,ispec,i,j,ier
-
- ! allocates temporary arrays
- allocate(mask_ibool(nglob),stat=ier)
- if( ier /= 0 ) stop 'error allocating mask_ibool'
- allocate(copy_ibool_ori(NGLLX,NGLLZ,nspec),stat=ier)
- if( ier /= 0 ) stop 'error allocating copy_ibool_ori'
+ integer :: inumber,ispec,i,j
! initializes temporary arrays
- mask_ibool(:) = -1
- copy_ibool_ori(:,:,:) = ibool(:,:,:)
+ integer_mask_ibool(:) = -1
inumber = 0
@@ -73,51 +67,42 @@
do ispec = 1,nspec
do j=1,NGLLZ
do i=1,NGLLX
- if(mask_ibool(copy_ibool_ori(i,j,ispec)) == -1) then
+ if(integer_mask_ibool(copy_ibool_ori(i,j,ispec)) == -1) then
! create a new point
inumber = inumber + 1
ibool(i,j,ispec) = inumber
- mask_ibool(copy_ibool_ori(i,j,ispec)) = inumber
+ integer_mask_ibool(copy_ibool_ori(i,j,ispec)) = inumber
else
! use an existing point created previously
- ibool(i,j,ispec) = mask_ibool(copy_ibool_ori(i,j,ispec))
+ ibool(i,j,ispec) = integer_mask_ibool(copy_ibool_ori(i,j,ispec))
endif
enddo
enddo
enddo
- deallocate(mask_ibool,copy_ibool_ori)
-
end subroutine get_global
-
!
!-------------------------------------------------------------------------------------------------
!
- subroutine get_global_indirect_addressing(nspec,nglob,ibool)
+! create a sorted version of the indirect addressing to reduce cache misses
-!- we can create a new indirect addressing to reduce cache misses
+ subroutine get_global_indirect_addressing(nspec,nglob,ibool,copy_ibool_ori,integer_mask_ibool)
implicit none
include "constants.h"
integer :: nspec,nglob
- integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
- ! local parameters
- integer, dimension(:,:,:), allocatable :: copy_ibool_ori
- integer, dimension(:), allocatable :: mask_ibool
- integer :: inumber,ispec,i,j,ier
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool,copy_ibool_ori
+ integer, dimension(nglob) :: integer_mask_ibool
- ! allocates temporary arrays
- allocate(mask_ibool(nglob),stat=ier)
- if( ier /= 0 ) stop 'error allocating mask_ibool'
- allocate(copy_ibool_ori(NGLLX,NGLLZ,nspec),stat=ier)
- if( ier /= 0 ) stop 'error allocating copy_ibool_ori'
+ ! local parameters
+ integer :: inumber,ispec,i,j
! initializes temporary arrays
- mask_ibool(:) = -1
+ integer_mask_ibool(:) = -1
copy_ibool_ori(:,:,:) = ibool(:,:,:)
inumber = 0
@@ -125,19 +110,18 @@
do ispec = 1,nspec
do j=1,NGLLZ
do i=1,NGLLX
- if(mask_ibool(copy_ibool_ori(i,j,ispec)) == -1) then
+ if(integer_mask_ibool(copy_ibool_ori(i,j,ispec)) == -1) then
! create a new point
inumber = inumber + 1
ibool(i,j,ispec) = inumber
- mask_ibool(copy_ibool_ori(i,j,ispec)) = inumber
+ integer_mask_ibool(copy_ibool_ori(i,j,ispec)) = inumber
else
! use an existing point created previously
- ibool(i,j,ispec) = mask_ibool(copy_ibool_ori(i,j,ispec))
+ ibool(i,j,ispec) = integer_mask_ibool(copy_ibool_ori(i,j,ispec))
endif
enddo
enddo
enddo
- deallocate(mask_ibool,copy_ibool_ori)
-
end subroutine get_global_indirect_addressing
+
diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90
index 4661082..fb10553 100644
--- a/src/specfem2D/specfem2D.F90
+++ b/src/specfem2D/specfem2D.F90
@@ -879,6 +879,10 @@
logical :: this_is_the_first_time_we_dump
logical, dimension(:), allocatable :: mask_ibool,mask_ibool_pml
+! to create a sorted version of the indirect addressing to reduce cache misses
+ integer, dimension(:,:,:), allocatable :: copy_ibool_ori
+ integer, dimension(:), allocatable :: integer_mask_ibool
+
double precision, dimension(:,:,:),allocatable:: rho_local,vp_local,vs_local
!!!! hessian
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rhorho_el_hessian_final1, rhorho_el_hessian_final2
@@ -896,6 +900,7 @@
! for horizontal periodic conditions
logical :: ADD_PERIODIC_CONDITIONS
+ logical, dimension(:), allocatable :: this_ibool_is_a_periodic_edge,copy_this_ibool_is_a_periodic
! horizontal periodicity distance for periodic conditions
double precision :: PERIODIC_HORIZ_DIST
@@ -1863,7 +1868,6 @@
x_center_spring = (xmax + xmin)/2.d0
z_center_spring = (zmax + zmin)/2.d0
-! YYYYYYYYYYYYYYYYYYYYYYYYyyyyyyyyyyyyyyyyyyyy
! periodic conditions: detect common points between left and right edges and replace one of them with the other
if(ADD_PERIODIC_CONDITIONS) then
@@ -1886,6 +1890,10 @@
stop 'periodic conditions currently implemented for a serial simulation only'
#endif
+! allocate an array to make sure that an acoustic free surface is not enforced on periodic edges
+ allocate(this_ibool_is_a_periodic_edge(NGLOB))
+ allocate(copy_this_ibool_is_a_periodic(NGLOB))
+
! set up a local geometric tolerance
xtypdist = +HUGEVAL
@@ -1931,6 +1939,7 @@
if (myrank == 0) write(IOUT,*) &
'start detecting points for periodic boundary conditions (the current algorithm can be slow and could be improved)...'
counter = 0
+ this_ibool_is_a_periodic_edge(:) = .false.
do iglob = 1,NGLOB-1
do iglob2 = iglob + 1,NGLOB
! check if the two points have the exact same Z coordinate
@@ -1940,6 +1949,8 @@
! if so, they are the same point, thus replace the highest value of ibool with the lowest
! to make them the same global point and thus implement periodicity automatically
counter = counter + 1
+ this_ibool_is_a_periodic_edge(iglob) = .true.
+ this_ibool_is_a_periodic_edge(iglob2) = .true.
do ispec = 1,nspec
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -1960,8 +1971,12 @@
if(counter > 0) write(IOUT,*) 'implemented periodic conditions on ',counter,' grid points on proc ',myrank
+ else
+
+ ! dummy allocation just to be able to use this array as a subroutine argument later
+ allocate(this_ibool_is_a_periodic_edge(1))
+
endif ! of if(ADD_PERIODIC_CONDITIONS)
-! end of periodic conditions
#ifdef USE_MPI
call MPI_REDUCE(count_nspec_acoustic, count_nspec_acoustic_total, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier)
@@ -2010,9 +2025,30 @@
endif
endif
+ ! allocate temporary arrays
+ allocate(integer_mask_ibool(nglob),stat=ier)
+ if( ier /= 0 ) stop 'error allocating mask_ibool'
+ allocate(copy_ibool_ori(NGLLX,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating copy_ibool_ori'
+
! reduce cache misses by sorting the global numbering in the order in which it is accessed in the time loop.
! this speeds up the calculations significantly on modern processors
- call get_global(nspec,nglob,ibool)
+ copy_ibool_ori(:,:,:) = ibool(:,:,:)
+ call get_global(nspec,nglob,ibool,copy_ibool_ori,integer_mask_ibool)
+
+ ! put the periodic edge flag in the new ibool order
+ if(ADD_PERIODIC_CONDITIONS) then
+ copy_this_ibool_is_a_periodic(:) = this_ibool_is_a_periodic_edge(:)
+ this_ibool_is_a_periodic_edge(:) = .false.
+ do ispec = 1,nspec
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ ! if this point is on a periodic edge in the old numbering system, set the flag in the new numbering system
+ if(copy_this_ibool_is_a_periodic(copy_ibool_ori(i,j,ispec))) this_ibool_is_a_periodic_edge(ibool(i,j,ispec)) = .true.
+ enddo
+ enddo
+ enddo
+ endif
!---- compute shape functions and their derivatives for regular interpolated display grid
do j = 1,pointsdisp
@@ -3239,13 +3275,13 @@
enddo
! reduces cache misses for outer elements
- call get_global_indirect_addressing(nspec_outer,nglob,ibool_outer)
+ call get_global_indirect_addressing(nspec_outer,nglob,ibool_outer,copy_ibool_ori,integer_mask_ibool)
! the total number of points without multiples in this region is now known
nglob_outer = maxval(ibool_outer)
! reduces cache misses for inner elements
- call get_global_indirect_addressing(nspec_inner,nglob,ibool_inner)
+ call get_global_indirect_addressing(nspec_inner,nglob,ibool_inner,copy_ibool_ori,integer_mask_ibool)
! the total number of points without multiples in this region is now known
nglob_inner = maxval(ibool_inner)
@@ -4714,12 +4750,12 @@ if(coupled_elastic_poro) then
if ( nelem_acoustic_surface > 0 ) then
call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,acoustic_surface, &
- ibool,nelem_acoustic_surface,nglob,nspec)
+ ibool,nelem_acoustic_surface,nglob,nspec,this_ibool_is_a_periodic_edge)
if(SIMULATION_TYPE == 3) then ! Adjoint calculation
call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
b_potential_acoustic,acoustic_surface, &
- ibool,nelem_acoustic_surface,nglob,nspec)
+ ibool,nelem_acoustic_surface,nglob,nspec,this_ibool_is_a_periodic_edge)
endif
endif
@@ -5216,12 +5252,12 @@ if(coupled_elastic_poro) then
if ( nelem_acoustic_surface > 0 ) then
call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,acoustic_surface, &
- ibool,nelem_acoustic_surface,nglob,nspec)
+ ibool,nelem_acoustic_surface,nglob,nspec,this_ibool_is_a_periodic_edge)
if(SIMULATION_TYPE == 3) then
call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
b_potential_acoustic,acoustic_surface, &
- ibool,nelem_acoustic_surface,nglob,nspec)
+ ibool,nelem_acoustic_surface,nglob,nspec,this_ibool_is_a_periodic_edge)
endif
endif
@@ -7367,7 +7403,7 @@ if(coupled_elastic_poro) then
if ( nelem_acoustic_surface > 0 ) then
call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
b_potential_acoustic,acoustic_surface, &
- ibool,nelem_acoustic_surface,nglob,nspec)
+ ibool,nelem_acoustic_surface,nglob,nspec,this_ibool_is_a_periodic_edge)
endif
endif
More information about the CIG-COMMITS
mailing list