[cig-commits] r17931 - seismo/2D/SPECFEM2D/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Feb 21 17:39:25 PST 2011
Author: danielpeter
Date: 2011-02-21 17:39:25 -0800 (Mon, 21 Feb 2011)
New Revision: 17931
Added:
seismo/2D/SPECFEM2D/trunk/get_MPI.F90
seismo/2D/SPECFEM2D/trunk/prepare_assemble_MPI.F90
seismo/2D/SPECFEM2D/trunk/sort_array_coordinates.F90
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile.in
seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
seismo/2D/SPECFEM2D/trunk/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
seismo/2D/SPECFEM2D/trunk/gmat01.f90
seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
seismo/2D/SPECFEM2D/trunk/prepare_color_image.F90
seismo/2D/SPECFEM2D/trunk/read_databases.f90
seismo/2D/SPECFEM2D/trunk/read_external_model.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
updates MPI routines; adds new files get_MPI.F90, prepare_assemble_MPI.F90, sort_array_coordinates.F90
Modified: seismo/2D/SPECFEM2D/trunk/Makefile.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile.in 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/Makefile.in 2011-02-22 01:39:25 UTC (rev 17931)
@@ -149,6 +149,7 @@
$O/define_shape_functions.o \
$O/enforce_acoustic_free_surface.o \
$O/exit_mpi.o \
+ $O/get_MPI.o \
$O/get_perm_cuthill_mckee.o \
$O/get_poroelastic_velocities.o \
$O/gll_library.o \
@@ -166,6 +167,7 @@
$O/plotgll.o \
$O/plotpost.o \
$O/prepare_absorb.o \
+ $O/prepare_assemble_MPI.o \
$O/prepare_color_image.o \
$O/prepare_initialfield.o \
$O/prepare_source_time_function.o \
@@ -175,6 +177,7 @@
$O/save_openDX_jacobian.o \
$O/set_sources.o \
$O/setup_sources_receivers.o \
+ $O/sort_array_coordinates.o \
$O/write_seismograms.o \
$O/specfem2D.o
@@ -474,3 +477,13 @@
$O/read_databases.o: read_databases.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/read_databases.o read_databases.f90
+
+$O/prepare_assemble_MPI.o: prepare_assemble_MPI.F90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/prepare_assemble_MPI.o prepare_assemble_MPI.F90
+
+$O/get_MPI.o: get_MPI.F90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/get_MPI.o get_MPI.F90
+
+$O/sort_array_coordinates.o: sort_array_coordinates.F90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/sort_array_coordinates.o sort_array_coordinates.F90
+
Modified: seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -48,285 +48,20 @@
! These subroutines are for the most part not used in the sequential version.
!
-!-----------------------------------------------
-! Determines the points that are on the interfaces with other partitions, to help
-! build the communication buffers, and determines which elements are considered 'inner'
-! (no points in common with other partitions) and 'outer' (at least one point in common
-! with neighbouring partitions).
-! We have both acoustic and (poro)elastic buffers, for coupling between acoustic and (poro)elastic elements
-! led us to have two sets of communications.
-!-----------------------------------------------
-subroutine prepare_assemble_MPI (nspec,ibool, &
- knods, ngnod, &
- npoin, elastic, poroelastic, &
- ninterface, max_interface_size, &
- my_nelmnts_neighbours, my_interfaces, &
- ibool_interfaces_acoustic, ibool_interfaces_elastic, ibool_interfaces_poroelastic, &
- nibool_interfaces_acoustic, nibool_interfaces_elastic, nibool_interfaces_poroelastic, &
- inum_interfaces_acoustic, inum_interfaces_elastic, inum_interfaces_poroelastic, &
- ninterface_acoustic, ninterface_elastic, ninterface_poroelastic, &
- mask_ispec_inner_outer &
- )
- implicit none
-
- include 'constants.h'
-
- integer, intent(in) :: nspec, npoin, ngnod
- logical, dimension(nspec), intent(in) :: elastic, poroelastic
- integer, dimension(ngnod,nspec), intent(in) :: knods
- integer, dimension(NGLLX,NGLLZ,nspec), intent(in) :: ibool
-
- integer :: ninterface
- integer :: max_interface_size
- integer, dimension(ninterface) :: my_nelmnts_neighbours
- integer, dimension(4,max_interface_size,ninterface) :: my_interfaces
- integer, dimension(NGLLX*max_interface_size,ninterface) :: &
- ibool_interfaces_acoustic,ibool_interfaces_elastic,ibool_interfaces_poroelastic
- integer, dimension(ninterface) :: &
- nibool_interfaces_acoustic,nibool_interfaces_elastic,nibool_interfaces_poroelastic
-
- integer, dimension(ninterface), intent(out) :: &
- inum_interfaces_acoustic, inum_interfaces_elastic, inum_interfaces_poroelastic
- integer, intent(out) :: ninterface_acoustic, ninterface_elastic, ninterface_poroelastic
-
- integer :: num_interface
- integer :: ispec_interface
-
- logical, dimension(nspec), intent(inout) :: mask_ispec_inner_outer
-
- logical, dimension(npoin) :: mask_ibool_acoustic
- logical, dimension(npoin) :: mask_ibool_elastic
- logical, dimension(npoin) :: mask_ibool_poroelastic
-
- integer :: ixmin, ixmax
- integer :: izmin, izmax
- integer, dimension(ngnod) :: n
- integer :: e1, e2
- integer :: type
- integer :: ispec
-
- integer :: k
- integer :: npoin_interface_acoustic
- integer :: npoin_interface_elastic
- integer :: npoin_interface_poroelastic
-
- integer :: ix,iz
-
- integer :: sens
-
- ibool_interfaces_acoustic(:,:) = 0
- nibool_interfaces_acoustic(:) = 0
- ibool_interfaces_elastic(:,:) = 0
- nibool_interfaces_elastic(:) = 0
- ibool_interfaces_poroelastic(:,:) = 0
- nibool_interfaces_poroelastic(:) = 0
-
- do num_interface = 1, ninterface
- npoin_interface_acoustic = 0
- npoin_interface_elastic = 0
- npoin_interface_poroelastic = 0
- mask_ibool_acoustic(:) = .false.
- mask_ibool_elastic(:) = .false.
- mask_ibool_poroelastic(:) = .false.
-
- do ispec_interface = 1, my_nelmnts_neighbours(num_interface)
- ispec = my_interfaces(1,ispec_interface,num_interface)
- type = my_interfaces(2,ispec_interface,num_interface)
- do k = 1, ngnod
- n(k) = knods(k,ispec)
- end do
- e1 = my_interfaces(3,ispec_interface,num_interface)
- e2 = my_interfaces(4,ispec_interface,num_interface)
-
- call get_edge(ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax, sens)
-
- do iz = izmin, izmax, sens
- do ix = ixmin, ixmax, sens
-
- if ( elastic(ispec) ) then
-
- if(.not. mask_ibool_elastic(ibool(ix,iz,ispec))) then
- mask_ibool_elastic(ibool(ix,iz,ispec)) = .true.
- npoin_interface_elastic = npoin_interface_elastic + 1
- ibool_interfaces_elastic(npoin_interface_elastic,num_interface)=&
- ibool(ix,iz,ispec)
- end if
- else if ( poroelastic(ispec) ) then
-
- if(.not. mask_ibool_poroelastic(ibool(ix,iz,ispec))) then
- mask_ibool_poroelastic(ibool(ix,iz,ispec)) = .true.
- npoin_interface_poroelastic = npoin_interface_poroelastic + 1
- ibool_interfaces_poroelastic(npoin_interface_poroelastic,num_interface)=&
- ibool(ix,iz,ispec)
- end if
- else
-
- if(.not. mask_ibool_acoustic(ibool(ix,iz,ispec))) then
- mask_ibool_acoustic(ibool(ix,iz,ispec)) = .true.
- npoin_interface_acoustic = npoin_interface_acoustic + 1
- ibool_interfaces_acoustic(npoin_interface_acoustic,num_interface)=&
- ibool(ix,iz,ispec)
- end if
- end if
- end do
- end do
-
- end do
- nibool_interfaces_acoustic(num_interface) = npoin_interface_acoustic
- nibool_interfaces_elastic(num_interface) = npoin_interface_elastic
- nibool_interfaces_poroelastic(num_interface) = npoin_interface_poroelastic
-
- do ispec = 1, nspec
- do iz = 1, NGLLZ
- do ix = 1, NGLLX
- if ( mask_ibool_acoustic(ibool(ix,iz,ispec)) &
- .or. mask_ibool_elastic(ibool(ix,iz,ispec)) &
- .or. mask_ibool_poroelastic(ibool(ix,iz,ispec)) ) then
- mask_ispec_inner_outer(ispec) = .true.
- endif
-
- enddo
- enddo
- enddo
-
- end do
-
- ninterface_acoustic = 0
- ninterface_elastic = 0
- ninterface_poroelastic = 0
- do num_interface = 1, ninterface
- if ( nibool_interfaces_acoustic(num_interface) > 0 ) then
- ninterface_acoustic = ninterface_acoustic + 1
- inum_interfaces_acoustic(ninterface_acoustic) = num_interface
- end if
- if ( nibool_interfaces_elastic(num_interface) > 0 ) then
- ninterface_elastic = ninterface_elastic + 1
- inum_interfaces_elastic(ninterface_elastic) = num_interface
- end if
- if ( nibool_interfaces_poroelastic(num_interface) > 0 ) then
- ninterface_poroelastic = ninterface_poroelastic + 1
- inum_interfaces_poroelastic(ninterface_poroelastic) = num_interface
- end if
- end do
-
-end subroutine prepare_assemble_MPI
-
-
-!-----------------------------------------------
-! Get the points (ixmin, ixmax, izmin and izmax) on an node/edge for one element.
-! 'sens' is used to have DO loops with increment equal to 'sens' (-/+1).
-!-----------------------------------------------
-subroutine get_edge ( ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax, sens )
-
- implicit none
-
- include "constants.h"
-
- integer, intent(in) :: ngnod
- integer, dimension(ngnod), intent(in) :: n
- integer, intent(in) :: type, e1, e2
- integer, intent(out) :: ixmin, ixmax, izmin, izmax
- integer, intent(out) :: sens
-
- if ( type == 1 ) then
- if ( e1 == n(1) ) then
- ixmin = 1
- ixmax = 1
- izmin = 1
- izmax = 1
- end if
- if ( e1 == n(2) ) then
- ixmin = NGLLX
- ixmax = NGLLX
- izmin = 1
- izmax = 1
- end if
- if ( e1 == n(3) ) then
- ixmin = NGLLX
- ixmax = NGLLX
- izmin = NGLLZ
- izmax = NGLLZ
- end if
- if ( e1 == n(4) ) then
- ixmin = 1
- ixmax = 1
- izmin = NGLLZ
- izmax = NGLLZ
- end if
- sens = 1
- else
- if ( e1 == n(1) ) then
- ixmin = 1
- izmin = 1
- if ( e2 == n(2) ) then
- ixmax = NGLLX
- izmax = 1
- sens = 1
- end if
- if ( e2 == n(4) ) then
- ixmax = 1
- izmax = NGLLZ
- sens = 1
- end if
- end if
- if ( e1 == n(2) ) then
- ixmin = NGLLX
- izmin = 1
- if ( e2 == n(3) ) then
- ixmax = NGLLX
- izmax = NGLLZ
- sens = 1
- end if
- if ( e2 == n(1) ) then
- ixmax = 1
- izmax = 1
- sens = -1
- end if
- end if
- if ( e1 == n(3) ) then
- ixmin = NGLLX
- izmin = NGLLZ
- if ( e2 == n(4) ) then
- ixmax = 1
- izmax = NGLLZ
- sens = -1
- end if
- if ( e2 == n(2) ) then
- ixmax = NGLLX
- izmax = 1
- sens = -1
- end if
- end if
- if ( e1 == n(4) ) then
- ixmin = 1
- izmin = NGLLZ
- if ( e2 == n(1) ) then
- ixmax = 1
- izmax = 1
- sens = -1
- end if
- if ( e2 == n(3) ) then
- ixmax = NGLLX
- izmax = NGLLZ
- sens = 1
- end if
- end if
- end if
-
-end subroutine get_edge
-
-
#ifdef USE_MPI
!-----------------------------------------------
! Assembling the mass matrix.
!-----------------------------------------------
-subroutine assemble_MPI_scalar(array_val1, array_val2, array_val3, array_val4,npoin, &
- ninterface, max_interface_size, max_ibool_interfaces_size_ac, max_ibool_interfaces_size_el, &
- max_ibool_interfaces_size_po, &
- ibool_interfaces_acoustic,ibool_interfaces_elastic, ibool_interfaces_poroelastic, &
- nibool_interfaces_acoustic,nibool_interfaces_elastic,nibool_interfaces_poroelastic,my_neighbours)
+ subroutine assemble_MPI_scalar(array_val1, array_val2, array_val3, array_val4,npoin, &
+ ninterface, max_interface_size, max_ibool_interfaces_size_ac, &
+ max_ibool_interfaces_size_el, &
+ max_ibool_interfaces_size_po, &
+ ibool_interfaces_acoustic,ibool_interfaces_elastic, &
+ ibool_interfaces_poroelastic, &
+ nibool_interfaces_acoustic,nibool_interfaces_elastic, &
+ nibool_interfaces_poroelastic,my_neighbours)
implicit none
@@ -336,11 +71,12 @@
integer, intent(in) :: npoin
integer, intent(in) :: ninterface
integer, intent(in) :: max_interface_size
- integer, intent(in) :: max_ibool_interfaces_size_ac,max_ibool_interfaces_size_el,max_ibool_interfaces_size_po
+ integer, intent(in) :: max_ibool_interfaces_size_ac,max_ibool_interfaces_size_el, &
+ max_ibool_interfaces_size_po
integer, dimension(NGLLX*max_interface_size,ninterface), intent(in) :: &
- ibool_interfaces_acoustic,ibool_interfaces_elastic,ibool_interfaces_poroelastic
+ ibool_interfaces_acoustic,ibool_interfaces_elastic,ibool_interfaces_poroelastic
integer, dimension(ninterface), intent(in) :: nibool_interfaces_acoustic,nibool_interfaces_elastic, &
- nibool_interfaces_poroelastic
+ nibool_interfaces_poroelastic
integer, dimension(ninterface), intent(in) :: my_neighbours
! array to assemble
real(kind=CUSTOM_REAL), dimension(npoin), intent(inout) :: array_val1,array_val2,array_val3,array_val4
@@ -355,6 +91,8 @@
integer :: msg_status(MPI_STATUS_SIZE)
integer, dimension(ninterface) :: msg_requests
+ buffer_send_faces_scalar(:,:) = 0.d0
+ buffer_recv_faces_scalar(:,:) = 0.d0
do num_interface = 1, ninterface
@@ -402,32 +140,36 @@
ipoin = 0
do i = 1, nibool_interfaces_acoustic(num_interface)
ipoin = ipoin + 1
- array_val1(ibool_interfaces_acoustic(i,num_interface)) = array_val1(ibool_interfaces_acoustic(i,num_interface)) + &
- buffer_recv_faces_scalar(ipoin,num_interface)
+ array_val1(ibool_interfaces_acoustic(i,num_interface)) = &
+ array_val1(ibool_interfaces_acoustic(i,num_interface)) &
+ + buffer_recv_faces_scalar(ipoin,num_interface)
end do
do i = 1, nibool_interfaces_elastic(num_interface)
ipoin = ipoin + 1
- array_val2(ibool_interfaces_elastic(i,num_interface)) = array_val2(ibool_interfaces_elastic(i,num_interface)) + &
- buffer_recv_faces_scalar(ipoin,num_interface)
+ array_val2(ibool_interfaces_elastic(i,num_interface)) = &
+ array_val2(ibool_interfaces_elastic(i,num_interface)) &
+ + buffer_recv_faces_scalar(ipoin,num_interface)
end do
do i = 1, nibool_interfaces_poroelastic(num_interface)
ipoin = ipoin + 1
- array_val3(ibool_interfaces_poroelastic(i,num_interface)) = array_val3(ibool_interfaces_poroelastic(i,num_interface)) + &
- buffer_recv_faces_scalar(ipoin,num_interface)
+ array_val3(ibool_interfaces_poroelastic(i,num_interface)) = &
+ array_val3(ibool_interfaces_poroelastic(i,num_interface)) &
+ + buffer_recv_faces_scalar(ipoin,num_interface)
end do
do i = 1, nibool_interfaces_poroelastic(num_interface)
ipoin = ipoin + 1
- array_val4(ibool_interfaces_poroelastic(i,num_interface)) = array_val4(ibool_interfaces_poroelastic(i,num_interface)) + &
- buffer_recv_faces_scalar(ipoin,num_interface)
+ array_val4(ibool_interfaces_poroelastic(i,num_interface)) = &
+ array_val4(ibool_interfaces_poroelastic(i,num_interface)) &
+ + buffer_recv_faces_scalar(ipoin,num_interface)
end do
end do
call MPI_BARRIER(mpi_comm_world,ier)
-end subroutine assemble_MPI_scalar
+ end subroutine assemble_MPI_scalar
!-----------------------------------------------
@@ -441,16 +183,15 @@
! Particular care should be taken concerning possible optimisations of the
! communication scheme.
!-----------------------------------------------
-subroutine assemble_MPI_vector_ac(array_val1,npoin, &
- ninterface, ninterface_acoustic, &
- inum_interfaces_acoustic, &
- max_interface_size, max_ibool_interfaces_size_ac,&
- ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
- tab_requests_send_recv_acoustic, &
- buffer_send_faces_vector_ac, &
- buffer_recv_faces_vector_ac, &
- my_neighbours &
- )
+ subroutine assemble_MPI_vector_ac(array_val1,npoin, &
+ ninterface, ninterface_acoustic, &
+ inum_interfaces_acoustic, &
+ max_interface_size, max_ibool_interfaces_size_ac,&
+ ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+ tab_requests_send_recv_acoustic, &
+ buffer_send_faces_vector_ac, &
+ buffer_recv_faces_vector_ac, &
+ my_neighbours )
implicit none
@@ -470,16 +211,20 @@
buffer_send_faces_vector_ac
real(kind=CUSTOM_REAL), dimension(max_ibool_interfaces_size_ac,ninterface_acoustic), intent(inout) :: &
buffer_recv_faces_vector_ac
- ! array to assemble
+ ! array to assemble
real(kind=CUSTOM_REAL), dimension(npoin), intent(inout) :: array_val1
integer, dimension(ninterface), intent(in) :: my_neighbours
+ ! local parameters
integer :: ipoin, num_interface, inum_interface
integer :: ier
integer :: i
integer, dimension(MPI_STATUS_SIZE) :: status_acoustic
-
+ ! initializes buffers
+ buffer_send_faces_vector_ac(:,:) = 0._CUSTOM_REAL
+ buffer_recv_faces_vector_ac(:,:) = 0._CUSTOM_REAL
+
do inum_interface = 1, ninterface_acoustic
num_interface = inum_interfaces_acoustic(inum_interface)
@@ -530,13 +275,14 @@
ipoin = 0
do i = 1, nibool_interfaces_acoustic(num_interface)
ipoin = ipoin + 1
- array_val1(ibool_interfaces_acoustic(i,num_interface)) = array_val1(ibool_interfaces_acoustic(i,num_interface)) + &
- buffer_recv_faces_vector_ac(ipoin,inum_interface)
+ array_val1(ibool_interfaces_acoustic(i,num_interface)) = &
+ array_val1(ibool_interfaces_acoustic(i,num_interface)) &
+ + buffer_recv_faces_vector_ac(ipoin,inum_interface)
end do
end do
-end subroutine assemble_MPI_vector_ac
+ end subroutine assemble_MPI_vector_ac
!-----------------------------------------------
@@ -550,7 +296,7 @@
! Particular care should be taken concerning possible optimisations of the
! communication scheme.
!-----------------------------------------------
-subroutine assemble_MPI_vector_el(array_val2,npoin, &
+ subroutine assemble_MPI_vector_el(array_val2,npoin, &
ninterface, ninterface_elastic, &
inum_interfaces_elastic, &
max_interface_size, max_ibool_interfaces_size_el,&
@@ -579,7 +325,7 @@
buffer_send_faces_vector_el
real(CUSTOM_REAL), dimension(max_ibool_interfaces_size_el,ninterface_elastic), intent(inout) :: &
buffer_recv_faces_vector_el
- ! array to assemble
+ ! array to assemble
real(kind=CUSTOM_REAL), dimension(3,npoin), intent(inout) :: array_val2
integer, dimension(ninterface), intent(in) :: my_neighbours
@@ -639,14 +385,15 @@
ipoin = 0
do i = 1, nibool_interfaces_elastic(num_interface)
- array_val2(:,ibool_interfaces_elastic(i,num_interface)) = array_val2(:,ibool_interfaces_elastic(i,num_interface)) + &
- buffer_recv_faces_vector_el(ipoin+1:ipoin+3,inum_interface)
+ array_val2(:,ibool_interfaces_elastic(i,num_interface)) = &
+ array_val2(:,ibool_interfaces_elastic(i,num_interface)) &
+ + buffer_recv_faces_vector_el(ipoin+1:ipoin+3,inum_interface)
ipoin = ipoin + 3
end do
end do
-end subroutine assemble_MPI_vector_el
+ end subroutine assemble_MPI_vector_el
!-----------------------------------------------
@@ -660,7 +407,7 @@
! Particular care should be taken concerning possible optimisations of the
! communication scheme.
!-----------------------------------------------
-subroutine assemble_MPI_vector_po(array_val3,array_val4,npoin, &
+ subroutine assemble_MPI_vector_po(array_val3,array_val4,npoin, &
ninterface, ninterface_poroelastic, &
inum_interfaces_poroelastic, &
max_interface_size, max_ibool_interfaces_size_po,&
@@ -689,7 +436,7 @@
buffer_send_faces_vector_pos,buffer_send_faces_vector_pow
real(CUSTOM_REAL), dimension(max_ibool_interfaces_size_po,ninterface_poroelastic), intent(inout) :: &
buffer_recv_faces_vector_pos,buffer_recv_faces_vector_pow
-! array to assemble
+ ! array to assemble
real(kind=CUSTOM_REAL), dimension(NDIM,npoin), intent(inout) :: array_val3,array_val4
integer, dimension(ninterface), intent(in) :: my_neighbours
@@ -790,6 +537,6 @@
end do
-end subroutine assemble_MPI_vector_po
+ end subroutine assemble_MPI_vector_po
#endif
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -44,8 +44,10 @@
subroutine compute_energy(displ_elastic,veloc_elastic,displs_poroelastic,velocs_poroelastic, &
displw_poroelastic,velocw_poroelastic, &
- xix,xiz,gammax,gammaz,jacobian,ibool,elastic,poroelastic,hprime_xx,hprime_zz, &
- nspec,npoin,assign_external_model,it,deltat,t0,kmato,elastcoef,density, &
+ xix,xiz,gammax,gammaz,jacobian,ibool, &
+ elastic,poroelastic,hprime_xx,hprime_zz, &
+ nspec,npoin,assign_external_model,it,deltat, &
+ t0,kmato,elastcoef,density, &
porosity,tortuosity, &
vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext, &
anisotropic,anisotropy,wxgll,wzgll,numat, &
@@ -71,7 +73,8 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
- real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_acoustic,potential_dot_dot_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin) :: &
+ potential_dot_acoustic,potential_dot_dot_acoustic
logical :: TURN_ATTENUATION_ON
@@ -93,7 +96,8 @@
double precision, dimension(6,numat) :: anisotropy
double precision, dimension(4,3,numat) :: elastcoef
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,vsext,rhoext
- double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext, &
+ c33ext,c35ext,c55ext
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displ_elastic,veloc_elastic
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,velocs_poroelastic
@@ -120,7 +124,8 @@
real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl,jacobianl
real(kind=CUSTOM_REAL) :: kinetic_energy,potential_energy
- real(kind=CUSTOM_REAL) :: cpl,csl,rhol,mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed,kappal
+ real(kind=CUSTOM_REAL) :: cpl,csl,rhol,mul_relaxed,lambdal_relaxed, &
+ lambdalplus2mul_relaxed,kappal
real(kind=CUSTOM_REAL) :: mul_s,kappal_s,rhol_s
real(kind=CUSTOM_REAL) :: kappal_f,rhol_f
real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr,phil,tortl
@@ -133,22 +138,22 @@
! loop over spectral elements
do ispec = 1,nspec
-!---
-!--- elastic spectral element
-!---
+ !---
+ !--- elastic spectral element
+ !---
if(elastic(ispec)) then
-! get relaxed elastic parameters of current spectral element
+ ! get relaxed elastic parameters of current spectral element
lambdal_relaxed = elastcoef(1,1,kmato(ispec))
mul_relaxed = elastcoef(2,1,kmato(ispec))
lambdalplus2mul_relaxed = elastcoef(3,1,kmato(ispec))
rhol = density(1,kmato(ispec))
-! double loop over GLL points
+ ! double loop over GLL points
do j = 1,NGLLZ
do i = 1,NGLLX
-!--- if external medium, get elastic parameters of current grid point
+ !--- if external medium, get elastic parameters of current grid point
if(assign_external_model) then
cpl = vpext(i,j,ispec)
csl = vsext(i,j,ispec)
@@ -158,15 +163,15 @@
lambdalplus2mul_relaxed = lambdal_relaxed + TWO*mul_relaxed
endif
-! derivative along x and along z
+ ! derivative along x and along z
dux_dxi = ZERO
duz_dxi = ZERO
dux_dgamma = ZERO
duz_dgamma = ZERO
-! first double loop over GLL points to compute and store gradients
-! we can merge the two loops because NGLLX == NGLLZ
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
duz_dxi = duz_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
@@ -180,64 +185,67 @@
gammazl = gammaz(i,j,ispec)
jacobianl = jacobian(i,j,ispec)
-! derivatives of displacement
+ ! derivatives of displacement
dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
-! compute kinetic energy
- kinetic_energy = kinetic_energy + &
- rhol*(veloc_elastic(1,ibool(i,j,ispec))**2 + &
- veloc_elastic(2,ibool(i,j,ispec))**2) *wxgll(i)*wzgll(j)*jacobianl / TWO
+ ! compute kinetic energy
+ kinetic_energy = kinetic_energy &
+ + rhol*(veloc_elastic(1,ibool(i,j,ispec))**2 &
+ + veloc_elastic(2,ibool(i,j,ispec))**2) *wxgll(i)*wzgll(j)*jacobianl / TWO
-! compute potential energy
- potential_energy = potential_energy + (lambdalplus2mul_relaxed*dux_dxl**2 &
+ ! compute potential energy
+ potential_energy = potential_energy &
+ + (lambdalplus2mul_relaxed*dux_dxl**2 &
+ lambdalplus2mul_relaxed*duz_dzl**2 &
- + two*lambdal_relaxed*dux_dxl*duz_dzl + mul_relaxed*(dux_dzl + duz_dxl)**2)*wxgll(i)*wzgll(j)*jacobianl / TWO
+ + two*lambdal_relaxed*dux_dxl*duz_dzl &
+ + mul_relaxed*(dux_dzl + duz_dxl)**2)*wxgll(i)*wzgll(j)*jacobianl / TWO
enddo
enddo
-!---
-!--- poroelastic spectral element
-!---
+ !---
+ !--- poroelastic spectral element
+ !---
elseif(poroelastic(ispec)) then
-! get relaxed elastic parameters of current spectral element
-!for now replaced by solid, fluid, and frame parameters of current spectral element
- phil = porosity(kmato(ispec))
- tortl = tortuosity(kmato(ispec))
-!solid properties
- mul_s = elastcoef(2,1,kmato(ispec))
- kappal_s = elastcoef(3,1,kmato(ispec)) - FOUR_THIRDS*mul_s
- rhol_s = density(1,kmato(ispec))
-!fluid properties
- kappal_f = elastcoef(1,2,kmato(ispec))
- rhol_f = density(2,kmato(ispec))
-!frame properties
- mul_fr = elastcoef(2,3,kmato(ispec))
- kappal_fr = elastcoef(3,3,kmato(ispec)) - FOUR_THIRDS*mul_fr
- rhol_bar = (1.d0 - phil)*rhol_s + phil*rhol_f
-!Biot coefficients for the input phi
+ ! get relaxed elastic parameters of current spectral element
+ !for now replaced by solid, fluid, and frame parameters of current spectral element
+ phil = porosity(kmato(ispec))
+ tortl = tortuosity(kmato(ispec))
+ !solid properties
+ mul_s = elastcoef(2,1,kmato(ispec))
+ kappal_s = elastcoef(3,1,kmato(ispec)) - FOUR_THIRDS*mul_s
+ rhol_s = density(1,kmato(ispec))
+ !fluid properties
+ kappal_f = elastcoef(1,2,kmato(ispec))
+ rhol_f = density(2,kmato(ispec))
+ !frame properties
+ mul_fr = elastcoef(2,3,kmato(ispec))
+ kappal_fr = elastcoef(3,3,kmato(ispec)) - FOUR_THIRDS*mul_fr
+ rhol_bar = (1.d0 - phil)*rhol_s + phil*rhol_f
+ !Biot coefficients for the input phi
D_biot = kappal_s*(1.d0 + phil*(kappal_s/kappal_f - 1.d0))
- H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) + kappal_fr + FOUR_THIRDS*mul_fr
+ H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) &
+ + kappal_fr + FOUR_THIRDS*mul_fr
C_biot = kappal_s*(kappal_s - kappal_fr)/(D_biot - kappal_fr)
M_biot = kappal_s*kappal_s/(D_biot - kappal_fr)
-!The RHS has the form : div T -phi/c div T_f + phi/ceta_fk^-1.partial t w
-!where T = G:grad u_s + C div w I
-!and T_f = C div u_s I + M div w I
-!we are expressing lambdaplus2mu, lambda, and mu for G, C, and M
+ !The RHS has the form : div T -phi/c div T_f + phi/ceta_fk^-1.partial t w
+ !where T = G:grad u_s + C div w I
+ !and T_f = C div u_s I + M div w I
+ !we are expressing lambdaplus2mu, lambda, and mu for G, C, and M
mul_G = mul_fr
lambdal_G = H_biot - TWO*mul_fr
lambdalplus2mul_G = lambdal_G + TWO*mul_G
-! first double loop over GLL points to compute and store gradients
+ ! first double loop over GLL points to compute and store gradients
do j = 1,NGLLZ
do i = 1,NGLLX
-! derivative along x and along z
+ ! derivative along x and along z
dux_dxi = ZERO
duz_dxi = ZERO
@@ -250,8 +258,8 @@
dwx_dgamma = ZERO
dwz_dgamma = ZERO
-! first double loop over GLL points to compute and store gradients
-! we can merge the two loops because NGLLX == NGLLZ
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
dux_dxi = dux_dxi + displs_poroelastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
duz_dxi = duz_dxi + displs_poroelastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
@@ -271,7 +279,7 @@
gammazl = gammaz(i,j,ispec)
jacobianl = jacobian(i,j,ispec)
-! derivatives of displacement
+ ! derivatives of displacement
dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
@@ -284,8 +292,9 @@
dwz_dxl = dwz_dxi*xixl + dwz_dgamma*gammaxl
dwz_dzl = dwz_dxi*xizl + dwz_dgamma*gammazl
-! compute potential energy
- potential_energy = potential_energy + ( lambdalplus2mul_G*dux_dxl**2 &
+ ! compute potential energy
+ potential_energy = potential_energy &
+ + ( lambdalplus2mul_G*dux_dxl**2 &
+ lambdalplus2mul_G*duz_dzl**2 &
+ two*lambdal_G*dux_dxl*duz_dzl + mul_G*(dux_dzl + duz_dxl)**2 &
+ two*C_biot*dwx_dxl*dux_dxl + two*C_biot*dwz_dzl*duz_dzl &
@@ -293,74 +302,84 @@
+ M_biot*dwx_dxl**2 + M_biot*dwz_dzl**2 &
+ two*M_biot*dwx_dxl*dwz_dzl )*wxgll(i)*wzgll(j)*jacobianl / TWO
-! compute kinetic energy
- if(phil > 0.0d0) then
- kinetic_energy = kinetic_energy + ( &
- rhol_bar*(velocs_poroelastic(1,ibool(i,j,ispec))**2 + velocs_poroelastic(2,ibool(i,j,ispec))**2) &
- + rhol_f*tortl/phil*(velocw_poroelastic(1,ibool(i,j,ispec))**2 + velocw_poroelastic(2,ibool(i,j,ispec))**2) &
+ ! compute kinetic energy
+ if(phil > 0.0d0) then
+ kinetic_energy = kinetic_energy &
+ + ( rhol_bar*(velocs_poroelastic(1,ibool(i,j,ispec))**2 &
+ + velocs_poroelastic(2,ibool(i,j,ispec))**2) &
+ + rhol_f*tortl/phil*(velocw_poroelastic(1,ibool(i,j,ispec))**2 &
+ + velocw_poroelastic(2,ibool(i,j,ispec))**2) &
+ rhol_f*(velocs_poroelastic(1,ibool(i,j,ispec))*velocw_poroelastic(1,ibool(i,j,ispec)) &
+ velocs_poroelastic(2,ibool(i,j,ispec))*velocw_poroelastic(2,ibool(i,j,ispec))) &
)*wxgll(i)*wzgll(j)*jacobianl / TWO
- else
- kinetic_energy = kinetic_energy + &
- rhol_s*(velocs_poroelastic(1,ibool(i,j,ispec))**2 + velocs_poroelastic(2,ibool(i,j,ispec))**2) &
- *wxgll(i)*wzgll(j)*jacobianl / TWO
- endif
+ else
+ kinetic_energy = kinetic_energy &
+ + rhol_s*(velocs_poroelastic(1,ibool(i,j,ispec))**2 &
+ + velocs_poroelastic(2,ibool(i,j,ispec))**2)*wxgll(i)*wzgll(j)*jacobianl / TWO
+ endif
enddo
enddo
-!---
-!--- acoustic spectral element
-!---
+ !---
+ !--- acoustic spectral element
+ !---
else
-! for the definition of potential energy in an acoustic fluid, see for instance
-! equation (23) of M. Maess et al., Journal of Sound and Vibration 296 (2006) 264-276
+ ! for the definition of potential energy in an acoustic fluid, see for instance
+ ! equation (23) of M. Maess et al., Journal of Sound and Vibration 296 (2006) 264-276
-! in case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette,
-! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential
-! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002).
-! This permits acoustic-elastic coupling based on a non-iterative time scheme.
-! Displacement is then: u = grad(Chi) / rho
-! Velocity is then: v = grad(Chi_dot) / rho (Chi_dot being the time derivative of Chi)
-! and pressure is: p = - Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi).
+ ! in case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette,
+ ! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential
+ ! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002).
+ ! This permits acoustic-elastic coupling based on a non-iterative time scheme.
+ ! Displacement is then: u = grad(Chi) / rho
+ ! Velocity is then: v = grad(Chi_dot) / rho (Chi_dot being the time derivative of Chi)
+ ! and pressure is: p = - Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi).
-! compute pressure in this element
- call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,elastcoef,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
- ispec,e1,e11, &
- TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
+ ! compute pressure in this element
+ call compute_pressure_one_element(pressure_element, &
+ potential_dot_dot_acoustic,displ_elastic,elastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,assign_external_model, &
+ numat,kmato,elastcoef,vpext,vsext,rhoext, &
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
+ ispec,e1,e11, &
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
-! compute velocity vector field in this element
- call compute_vector_one_element(vector_field_element,potential_dot_acoustic,veloc_elastic,elastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
+ ! compute velocity vector field in this element
+ call compute_vector_one_element(vector_field_element, &
+ potential_dot_acoustic,veloc_elastic,elastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
-! get density of current spectral element
+ ! get density of current spectral element
lambdal_relaxed = elastcoef(1,1,kmato(ispec))
mul_relaxed = elastcoef(2,1,kmato(ispec))
rhol = density(1,kmato(ispec))
kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
cpl = sqrt((kappal + 4._CUSTOM_REAL*mul_relaxed/3._CUSTOM_REAL)/rhol)
-! double loop over GLL points
+ ! double loop over GLL points
do j = 1,NGLLZ
do i = 1,NGLLX
-!--- if external medium, get density of current grid point
+ !--- if external medium, get density of current grid point
if(assign_external_model) then
cpl = vpext(i,j,ispec)
rhol = rhoext(i,j,ispec)
endif
-! compute kinetic energy
- kinetic_energy = kinetic_energy + &
- rhol*(vector_field_element(1,i,j)**2 + &
- vector_field_element(2,i,j)**2) *wxgll(i)*wzgll(j)*jacobianl / TWO
+ jacobianl = jacobian(i,j,ispec)
-! compute potential energy
- potential_energy = potential_energy + (pressure_element(i,j)**2)*wxgll(i)*wzgll(j)*jacobianl / (TWO * rhol * cpl**2)
+ ! compute kinetic energy
+ kinetic_energy = kinetic_energy &
+ + rhol*(vector_field_element(1,i,j)**2 &
+ + vector_field_element(2,i,j)**2) *wxgll(i)*wzgll(j)*jacobianl / TWO
+ ! compute potential energy
+ potential_energy = potential_energy &
+ + (pressure_element(i,j)**2)*wxgll(i)*wzgll(j)*jacobianl / (TWO * rhol * cpl**2)
+
enddo
enddo
@@ -368,7 +387,7 @@
enddo
-! save kinetic, potential and total energy for this time step in external file
+ ! save kinetic, potential and total energy for this time step in external file
write(IOUT_ENERGY,*) real(dble(it-1)*deltat - t0,4),real(kinetic_energy,4), &
real(potential_energy,4),real(kinetic_energy + potential_energy,4)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -157,8 +157,8 @@
dux_dgamma = dux_dgamma + potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
if(SIMULATION_TYPE == 2) then
- b_dux_dxi = b_dux_dxi + b_potential_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
- b_dux_dgamma = b_dux_dgamma + b_potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
+ b_dux_dxi = b_dux_dxi + b_potential_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
+ b_dux_dgamma = b_dux_dgamma + b_potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
endif
enddo
@@ -172,24 +172,24 @@
dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
if(SIMULATION_TYPE == 2) then
- b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
- b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
+ b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
+ b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
endif
jacobianl = jacobian(i,j,ispec)
! if external density model
- if(assign_external_model) rhol = rhoext(i,j,ispec)
+ if(assign_external_model) rhol = rhoext(i,j,ispec)
! for acoustic medium
! also add GLL integration weights
tempx1(i,j) = wzgll(j)*jacobianl*(xixl*dux_dxl + xizl*dux_dzl) / rhol
tempx2(i,j) = wxgll(i)*jacobianl*(gammaxl*dux_dxl + gammazl*dux_dzl) / rhol
- if(SIMULATION_TYPE == 2) then
- b_tempx1(i,j) = wzgll(j)*jacobianl*(xixl*b_dux_dxl + xizl*b_dux_dzl) /rhol
- b_tempx2(i,j) = wxgll(i)*jacobianl*(gammaxl*b_dux_dxl + gammazl*b_dux_dzl) /rhol
- endif
+ if(SIMULATION_TYPE == 2) then
+ b_tempx1(i,j) = wzgll(j)*jacobianl*(xixl*b_dux_dxl + xizl*b_dux_dzl) /rhol
+ b_tempx2(i,j) = wxgll(i)*jacobianl*(gammaxl*b_dux_dxl + gammazl*b_dux_dzl) /rhol
+ endif
enddo
enddo
@@ -264,13 +264,16 @@
! Sommerfeld condition if acoustic
if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
- b_absorb_acoustic_left(j,ib_left(ispecabs),it) = potential_dot_acoustic(iglob)*weight/cpl/rhol
+ b_absorb_acoustic_left(j,ib_left(ispecabs),it) = &
+ potential_dot_acoustic(iglob)*weight/cpl/rhol
elseif(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) - &
- b_absorb_acoustic_left(j,ib_left(ispecabs),NSTEP-it+1)
+ b_potential_dot_dot_acoustic(iglob) = &
+ b_potential_dot_dot_acoustic(iglob) - &
+ b_absorb_acoustic_left(j,ib_left(ispecabs),NSTEP-it+1)
endif
endif
@@ -304,14 +307,17 @@
! Sommerfeld condition if acoustic
if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
- b_absorb_acoustic_right(j,ib_right(ispecabs),it) = potential_dot_acoustic(iglob)*weight/cpl/rhol
+ b_absorb_acoustic_right(j,ib_right(ispecabs),it) = &
+ potential_dot_acoustic(iglob)*weight/cpl/rhol
elseif(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) - &
- b_absorb_acoustic_right(j,ib_right(ispecabs),NSTEP-it+1)
+ b_potential_dot_dot_acoustic(iglob) = &
+ b_potential_dot_dot_acoustic(iglob) - &
+ b_absorb_acoustic_right(j,ib_right(ispecabs),NSTEP-it+1)
endif
endif
@@ -349,13 +355,16 @@
! Sommerfeld condition if acoustic
if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
- b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),it) = potential_dot_acoustic(iglob)*weight/cpl/rhol
+ b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),it) = &
+ potential_dot_acoustic(iglob)*weight/cpl/rhol
elseif(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) - &
- b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),NSTEP-it+1)
+ b_potential_dot_dot_acoustic(iglob) = &
+ b_potential_dot_dot_acoustic(iglob) - &
+ b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),NSTEP-it+1)
endif
endif
@@ -393,13 +402,16 @@
! Sommerfeld condition if acoustic
if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
- b_absorb_acoustic_top(i,ib_top(ispecabs),it) = potential_dot_acoustic(iglob)*weight/cpl/rhol
+ b_absorb_acoustic_top(i,ib_top(ispecabs),it) = &
+ potential_dot_acoustic(iglob)*weight/cpl/rhol
elseif(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) - &
- b_absorb_acoustic_top(i,ib_top(ispecabs),NSTEP-it+1)
+ b_potential_dot_dot_acoustic(iglob) = &
+ b_potential_dot_dot_acoustic(iglob) - &
+ b_absorb_acoustic_top(i,ib_top(ispecabs),NSTEP-it+1)
endif
endif
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -218,7 +218,7 @@
if(elastic(ispec)) then
-! get relaxed elastic parameters of current spectral element
+ ! get relaxed elastic parameters of current spectral element
lambdal_relaxed = elastcoef(1,1,kmato(ispec))
mul_relaxed = elastcoef(2,1,kmato(ispec))
lambdalplus2mul_relaxed = elastcoef(3,1,kmato(ispec))
@@ -226,24 +226,24 @@
do j = 1,NGLLZ
do i = 1,NGLLX
-!--- if external medium, get elastic parameters of current grid point
+ !--- if external medium, get elastic parameters of current grid point
if(assign_external_model) then
cpl = vpext(i,j,ispec)
csl = vsext(i,j,ispec)
denst = rhoext(i,j,ispec)
mul_relaxed = denst*csl*csl
lambdal_relaxed = denst*cpl*cpl - TWO*mul_relaxed
- endif
+ endif
-! derivative along x and along z
+ ! derivative along x and along z
dux_dxi = ZERO
duz_dxi = ZERO
dux_dgamma = ZERO
duz_dgamma = ZERO
-! first double loop over GLL points to compute and store gradients
-! we can merge the two loops because NGLLX == NGLLZ
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
duz_dxi = duz_dxi + displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
@@ -256,13 +256,13 @@
gammaxl = gammax(i,j,ispec)
gammazl = gammaz(i,j,ispec)
-! derivatives of displacement
+ ! derivatives of displacement
dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
! compute diagonal components of the stress tensor (include attenuation or anisotropy if needed)
- if(TURN_ATTENUATION_ON) then
+ if(TURN_ATTENUATION_ON) then
! attenuation is implemented following the memory variable formulation of
! J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
@@ -270,61 +270,67 @@
! J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation in a linear
! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611 (1988).
-! compute unrelaxed elastic coefficients from formulas in Carcione 1993 page 111
- lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1(i,j,ispec) - mul_relaxed * Mu_nu2(i,j,ispec)
- mul_unrelaxed = mul_relaxed * Mu_nu2(i,j,ispec)
- lambdalplus2mul_unrelaxed = lambdal_unrelaxed + TWO*mul_unrelaxed
+ ! compute unrelaxed elastic coefficients from formulas in Carcione 1993 page 111
+ lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1(i,j,ispec) &
+ - mul_relaxed * Mu_nu2(i,j,ispec)
+ mul_unrelaxed = mul_relaxed * Mu_nu2(i,j,ispec)
+ lambdalplus2mul_unrelaxed = lambdal_unrelaxed + TWO*mul_unrelaxed
-! compute the stress using the unrelaxed Lame parameters (Carcione 1993, page 111)
- sigma_xx = lambdalplus2mul_unrelaxed*dux_dxl + lambdal_unrelaxed*duz_dzl
- sigma_zz = lambdalplus2mul_unrelaxed*duz_dzl + lambdal_unrelaxed*dux_dxl
+ ! compute the stress using the unrelaxed Lame parameters (Carcione 1993, page 111)
+ sigma_xx = lambdalplus2mul_unrelaxed*dux_dxl + lambdal_unrelaxed*duz_dzl
+ sigma_zz = lambdalplus2mul_unrelaxed*duz_dzl + lambdal_unrelaxed*dux_dxl
-! add the memory variables using the relaxed parameters (Carcione 1993, page 111)
-! beware: there is a bug in Carcione's equation (2c) for sigma_zz, we fixed it in the code below
- e1_sum = 0._CUSTOM_REAL
- e11_sum = 0._CUSTOM_REAL
+ ! add the memory variables using the relaxed parameters (Carcione 1993, page 111)
+ ! beware: there is a bug in Carcione's equation (2c) for sigma_zz, we fixed it in the code below
+ e1_sum = 0._CUSTOM_REAL
+ e11_sum = 0._CUSTOM_REAL
- do i_sls = 1,N_SLS
- e1_sum = e1_sum + e1(i,j,ispec,i_sls)
- e11_sum = e11_sum + e11(i,j,ispec,i_sls)
- enddo
+ do i_sls = 1,N_SLS
+ e1_sum = e1_sum + e1(i,j,ispec,i_sls)
+ e11_sum = e11_sum + e11(i,j,ispec,i_sls)
+ enddo
- sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed) * e1_sum + TWO * mul_relaxed * e11_sum
- sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed) * e1_sum - TWO * mul_relaxed * e11_sum
+ sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed) * e1_sum &
+ + TWO * mul_relaxed * e11_sum
+ sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed) * e1_sum &
+ - TWO * mul_relaxed * e11_sum
- else
+ else
-! no attenuation
- sigma_xx = lambdalplus2mul_relaxed*dux_dxl + lambdal_relaxed*duz_dzl
- sigma_zz = lambdalplus2mul_relaxed*duz_dzl + lambdal_relaxed*dux_dxl
+ ! no attenuation
+ sigma_xx = lambdalplus2mul_relaxed*dux_dxl + lambdal_relaxed*duz_dzl
+ sigma_zz = lambdalplus2mul_relaxed*duz_dzl + lambdal_relaxed*dux_dxl
- endif
+ endif
-! full anisotropy
- if(anisotropic(ispec)) then
- if(assign_external_model) then
- c11 = c11ext(i,j,ispec)
- c15 = c15ext(i,j,ispec)
- c13 = c13ext(i,j,ispec)
- c33 = c33ext(i,j,ispec)
- c35 = c35ext(i,j,ispec)
- c55 = c55ext(i,j,ispec)
- else
- c11 = anisotropy(1,kmato(ispec))
- c13 = anisotropy(2,kmato(ispec))
- c15 = anisotropy(3,kmato(ispec))
- c33 = anisotropy(4,kmato(ispec))
- c35 = anisotropy(5,kmato(ispec))
- c55 = anisotropy(6,kmato(ispec))
- endif
+ ! full anisotropy
+ if(anisotropic(ispec)) then
+ if(assign_external_model) then
+ c11 = c11ext(i,j,ispec)
+ c15 = c15ext(i,j,ispec)
+ c13 = c13ext(i,j,ispec)
+ c33 = c33ext(i,j,ispec)
+ c35 = c35ext(i,j,ispec)
+ c55 = c55ext(i,j,ispec)
+ else
+ c11 = anisotropy(1,kmato(ispec))
+ c13 = anisotropy(2,kmato(ispec))
+ c15 = anisotropy(3,kmato(ispec))
+ c33 = anisotropy(4,kmato(ispec))
+ c35 = anisotropy(5,kmato(ispec))
+ c55 = anisotropy(6,kmato(ispec))
+ endif
-! implement anisotropy in 2D
- sigma_xx = c11*dux_dxl + c15*(duz_dxl + dux_dzl) + c13*duz_dzl
- sigma_zz = c13*dux_dxl + c35*(duz_dxl + dux_dzl) + c33*duz_dzl
+ duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
+ dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
- endif
+ ! implement anisotropy in 2D
+ sigma_xx = c11*dux_dxl + c15*(duz_dxl + dux_dzl) + c13*duz_dzl
+ sigma_zz = c13*dux_dxl + c35*(duz_dxl + dux_dzl) + c33*duz_dzl
-! store pressure
+ endif
+
+ ! store pressure
pressure_element(i,j) = - (sigma_xx + sigma_zz) / 2.d0
enddo
@@ -332,36 +338,40 @@
elseif(poroelastic(ispec)) then
-! get poroelastic parameters of current spectral element
+ lambdal_relaxed = elastcoef(1,1,kmato(ispec))
+ mul_relaxed = elastcoef(2,1,kmato(ispec))
+
+ ! get poroelastic parameters of current spectral element
phil = porosity(kmato(ispec))
tortl = tortuosity(kmato(ispec))
-!solid properties
+ !solid properties
mul_s = elastcoef(2,1,kmato(ispec))
kappal_s = elastcoef(3,1,kmato(ispec)) - FOUR_THIRDS*mul_s
rhol_s = density(1,kmato(ispec))
-!fluid properties
+ !fluid properties
kappal_f = elastcoef(1,2,kmato(ispec))
rhol_f = density(2,kmato(ispec))
-!frame properties
+ !frame properties
mul_fr = elastcoef(2,3,kmato(ispec))
kappal_fr = elastcoef(3,3,kmato(ispec)) - FOUR_THIRDS*mul_fr
rhol_bar = (1.d0 - phil)*rhol_s + phil*rhol_f
-!Biot coefficients for the input phi
- D_biot = kappal_s*(1.d0 + phil*(kappal_s/kappal_f - 1.d0))
- H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) + kappal_fr + FOUR_THIRDS*mul_fr
- C_biot = kappal_s*(kappal_s - kappal_fr)/(D_biot - kappal_fr)
- M_biot = kappal_s*kappal_s/(D_biot - kappal_fr)
-!where T = G:grad u_s + C div w I
-!and T_f = C div u_s I + M div w I
-!we are expressing lambdaplus2mu, lambda, and mu for G, C, and M
- mul_G = mul_fr
- lambdal_G = H_biot - TWO*mul_fr
- lambdalplus2mul_G = lambdal_G + TWO*mul_G
+ !Biot coefficients for the input phi
+ D_biot = kappal_s*(1.d0 + phil*(kappal_s/kappal_f - 1.d0))
+ H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) &
+ + kappal_fr + FOUR_THIRDS*mul_fr
+ C_biot = kappal_s*(kappal_s - kappal_fr)/(D_biot - kappal_fr)
+ M_biot = kappal_s*kappal_s/(D_biot - kappal_fr)
+ !where T = G:grad u_s + C div w I
+ !and T_f = C div u_s I + M div w I
+ !we are expressing lambdaplus2mu, lambda, and mu for G, C, and M
+ mul_G = mul_fr
+ lambdal_G = H_biot - TWO*mul_fr
+ lambdalplus2mul_G = lambdal_G + TWO*mul_G
do j = 1,NGLLZ
do i = 1,NGLLX
-! derivative along x and along z
+ ! derivative along x and along z
dux_dxi = ZERO
duz_dxi = ZERO
@@ -374,8 +384,8 @@
dwx_dgamma = ZERO
dwz_dgamma = ZERO
-! first double loop over GLL points to compute and store gradients
-! we can merge the two loops because NGLLX == NGLLZ
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the two loops because NGLLX == NGLLZ
do k = 1,NGLLX
dux_dxi = dux_dxi + displs_poroelastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
duz_dxi = duz_dxi + displs_poroelastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
@@ -394,7 +404,7 @@
gammaxl = gammax(i,j,ispec)
gammazl = gammaz(i,j,ispec)
-! derivatives of displacement
+ ! derivatives of displacement
dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
@@ -403,7 +413,7 @@
! compute diagonal components of the stress tensor (include attenuation if needed)
- if(TURN_ATTENUATION_ON) then
+ if(TURN_ATTENUATION_ON) then
!-------------------- ATTENTION TO BE DEFINED ------------------------------!
! attenuation is implemented following the memory variable formulation of
@@ -412,39 +422,42 @@
! J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation in a linear
! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611 (1988).
-! compute unrelaxed elastic coefficients from formulas in Carcione 1993 page 111
- lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1(i,j,ispec) - mul_relaxed * Mu_nu2(i,j,ispec)
- mul_unrelaxed = mul_relaxed * Mu_nu2(i,j,ispec)
- lambdalplus2mul_unrelaxed = lambdal_unrelaxed + TWO*mul_unrelaxed
+ ! compute unrelaxed elastic coefficients from formulas in Carcione 1993 page 111
+ lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1(i,j,ispec) &
+ - mul_relaxed * Mu_nu2(i,j,ispec)
+ mul_unrelaxed = mul_relaxed * Mu_nu2(i,j,ispec)
+ lambdalplus2mul_unrelaxed = lambdal_unrelaxed + TWO*mul_unrelaxed
-! compute the stress using the unrelaxed Lame parameters (Carcione 1993, page 111)
- sigma_xx = lambdalplus2mul_unrelaxed*dux_dxl + lambdal_unrelaxed*duz_dzl
- sigma_zz = lambdalplus2mul_unrelaxed*duz_dzl + lambdal_unrelaxed*dux_dxl
+ ! compute the stress using the unrelaxed Lame parameters (Carcione 1993, page 111)
+ sigma_xx = lambdalplus2mul_unrelaxed*dux_dxl + lambdal_unrelaxed*duz_dzl
+ sigma_zz = lambdalplus2mul_unrelaxed*duz_dzl + lambdal_unrelaxed*dux_dxl
-! add the memory variables using the relaxed parameters (Carcione 1993, page 111)
-! beware: there is a bug in Carcione's equation (2c) for sigma_zz, we fixed it in the code below
- e1_sum = 0._CUSTOM_REAL
- e11_sum = 0._CUSTOM_REAL
+ ! add the memory variables using the relaxed parameters (Carcione 1993, page 111)
+ ! beware: there is a bug in Carcione's equation (2c) for sigma_zz, we fixed it in the code below
+ e1_sum = 0._CUSTOM_REAL
+ e11_sum = 0._CUSTOM_REAL
- do i_sls = 1,N_SLS
- e1_sum = e1_sum + e1(i,j,ispec,i_sls)
- e11_sum = e11_sum + e11(i,j,ispec,i_sls)
- enddo
+ do i_sls = 1,N_SLS
+ e1_sum = e1_sum + e1(i,j,ispec,i_sls)
+ e11_sum = e11_sum + e11(i,j,ispec,i_sls)
+ enddo
- sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed) * e1_sum + TWO * mul_relaxed * e11_sum
- sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed) * e1_sum - TWO * mul_relaxed * e11_sum
+ sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed) * e1_sum &
+ + TWO * mul_relaxed * e11_sum
+ sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed) * e1_sum &
+ - TWO * mul_relaxed * e11_sum
- else
+ else
-! no attenuation
- sigma_xx = lambdalplus2mul_G*dux_dxl + lambdal_G*duz_dzl + C_biot*(dwx_dxl + dwz_dzl)
- sigma_zz = lambdalplus2mul_G*duz_dzl + lambdal_G*dux_dxl + C_biot*(dwx_dxl + dwz_dzl)
+ ! no attenuation
+ sigma_xx = lambdalplus2mul_G*dux_dxl + lambdal_G*duz_dzl + C_biot*(dwx_dxl + dwz_dzl)
+ sigma_zz = lambdalplus2mul_G*duz_dzl + lambdal_G*dux_dxl + C_biot*(dwx_dxl + dwz_dzl)
- sigmap = C_biot*(dux_dxl + duz_dzl) + M_biot*(dwx_dxl + dwz_dzl)
+ sigmap = C_biot*(dux_dxl + duz_dzl) + M_biot*(dwx_dxl + dwz_dzl)
- endif
+ endif
-! store pressure
+ ! store pressure
pressure_element(i,j) = - (sigma_xx + sigma_zz) / 2.d0
! pressure_element2(i,j) = - sigmap
enddo
@@ -458,7 +471,7 @@
iglob = ibool(i,j,ispec)
-! store pressure
+ ! store pressure
pressure_element(i,j) = - potential_dot_dot_acoustic(iglob)
enddo
Added: seismo/2D/SPECFEM2D/trunk/get_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/get_MPI.F90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/get_MPI.F90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -0,0 +1,319 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 6.1
+! ------------------------------
+!
+! Copyright Universite de Pau, CNRS and INRIA, France,
+! and Princeton University / California Institute of Technology, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+#ifdef USE_MPI
+
+ subroutine get_MPI(nspec,ibool,knods,ngnod,npoin,elastic,poroelastic, &
+ ninterface, max_interface_size, &
+ my_nelmnts_neighbours,my_interfaces,my_neighbours, &
+ ibool_interfaces_acoustic, ibool_interfaces_elastic, &
+ ibool_interfaces_poroelastic, &
+ nibool_interfaces_acoustic, nibool_interfaces_elastic, &
+ nibool_interfaces_poroelastic, &
+ inum_interfaces_acoustic, inum_interfaces_elastic, &
+ inum_interfaces_poroelastic, &
+ ninterface_acoustic, ninterface_elastic, ninterface_poroelastic, &
+ mask_ispec_inner_outer, &
+ myrank,ipass,coord)
+
+! sets up the MPI interface for communication between partitions
+
+ implicit none
+
+ include "constants.h"
+ include 'mpif.h'
+
+ integer, intent(in) :: nspec, npoin, ngnod
+ logical, dimension(nspec), intent(in) :: elastic, poroelastic
+ integer, dimension(ngnod,nspec), intent(in) :: knods
+ integer, dimension(NGLLX,NGLLZ,nspec), intent(in) :: ibool
+
+ integer :: ninterface
+ integer :: max_interface_size
+ integer, dimension(ninterface) :: my_nelmnts_neighbours,my_neighbours
+ integer, dimension(4,max_interface_size,ninterface) :: my_interfaces
+
+ integer, dimension(NGLLX*max_interface_size,ninterface) :: &
+ ibool_interfaces_acoustic,ibool_interfaces_elastic,ibool_interfaces_poroelastic
+ integer, dimension(ninterface) :: &
+ nibool_interfaces_acoustic,nibool_interfaces_elastic,nibool_interfaces_poroelastic
+ integer, dimension(ninterface), intent(out) :: &
+ inum_interfaces_acoustic, inum_interfaces_elastic, inum_interfaces_poroelastic
+ integer, intent(out) :: ninterface_acoustic, ninterface_elastic, ninterface_poroelastic
+
+ logical, dimension(nspec), intent(inout) :: mask_ispec_inner_outer
+
+ integer :: myrank,ipass
+ double precision, dimension(NDIM,npoin) :: coord
+
+ !local parameters
+ double precision, dimension(:), allocatable :: xp,zp
+ double precision, dimension(:), allocatable :: work
+ integer, dimension(:), allocatable :: locval
+ integer, dimension(:), allocatable :: nibool_interfaces_true
+ ! for MPI buffers
+ integer, dimension(:), allocatable :: reorder_interface,ind,ninseg,iwork
+ integer, dimension(:), allocatable :: ibool_dummy
+! integer, dimension(:,:), allocatable :: ibool_interfaces_dummy
+ logical, dimension(:), allocatable :: ifseg
+ integer :: iinterface,ilocnum
+ integer :: num_points1, num_points2
+ ! assembly test
+ integer :: i,j,ispec,iglob,count,inum,ier,idomain
+ integer :: max_nibool_interfaces,num_nibool,num_interface
+ real(kind=CUSTOM_REAL), dimension(:),allocatable :: test_flag_cr
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: buffer_send_faces_vector_ac
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: buffer_recv_faces_vector_ac
+ integer, dimension(:), allocatable :: tab_requests_send_recv_acoustic
+
+ ! gets global indices for points on MPI interfaces
+ ! (defined by my_interfaces) between different partitions
+ ! and stores them in ibool_interfaces*** & nibool_interfaces*** (number of total points)
+ call prepare_assemble_MPI(nspec,ibool,knods, ngnod,npoin, elastic, poroelastic, &
+ ninterface, max_interface_size, &
+ my_nelmnts_neighbours, my_interfaces, &
+ ibool_interfaces_acoustic, ibool_interfaces_elastic, &
+ ibool_interfaces_poroelastic, &
+ nibool_interfaces_acoustic, nibool_interfaces_elastic, &
+ nibool_interfaces_poroelastic, &
+ inum_interfaces_acoustic, inum_interfaces_elastic, &
+ inum_interfaces_poroelastic, &
+ ninterface_acoustic, ninterface_elastic, ninterface_poroelastic, &
+ mask_ispec_inner_outer )
+
+
+ ! sorts ibool comm buffers lexicographically for all MPI interfaces
+ num_points1 = 0
+ num_points2 = 0
+ allocate(nibool_interfaces_true(ninterface))
+
+ do idomain = 1,3
+
+ ! checks number of interface in this domain
+ num_interface = 0
+ if( idomain == 1 ) then
+ num_interface = ninterface_acoustic
+ elseif( idomain == 2 ) then
+ num_interface = ninterface_elastic
+ elseif( idomain == 3 ) then
+ num_interface = ninterface_poroelastic
+ endif
+ if( num_interface == 0 ) cycle
+
+ ! loops over interfaces
+ do iinterface = 1, ninterface
+
+ ! number of global points in this interface
+ num_nibool = 0
+ if( idomain == 1 ) then
+ num_nibool = nibool_interfaces_acoustic(iinterface)
+ elseif( idomain == 2 ) then
+ num_nibool = nibool_interfaces_elastic(iinterface)
+ elseif( idomain == 3 ) then
+ num_nibool = nibool_interfaces_poroelastic(iinterface)
+ endif
+ ! checks if anything to sort
+ if( num_nibool == 0 ) cycle
+
+ allocate(xp(num_nibool))
+ allocate(zp(num_nibool))
+ allocate(locval(num_nibool))
+ allocate(ifseg(num_nibool))
+ allocate(reorder_interface(num_nibool))
+ allocate(ibool_dummy(num_nibool))
+ allocate(ind(num_nibool))
+ allocate(ninseg(num_nibool))
+ allocate(iwork(num_nibool))
+ allocate(work(num_nibool))
+
+ ! works with a copy of ibool array
+ if( idomain == 1 ) then
+ ibool_dummy(:) = ibool_interfaces_acoustic(1:num_nibool,iinterface)
+ elseif( idomain == 2 ) then
+ ibool_dummy(:) = ibool_interfaces_elastic(1:num_nibool,iinterface)
+ elseif( idomain == 3 ) then
+ ibool_dummy(:) = ibool_interfaces_poroelastic(1:num_nibool,iinterface)
+ endif
+
+ ! gets x,y,z coordinates of global points on MPI interface
+ do ilocnum = 1, num_nibool
+ iglob = ibool_dummy(ilocnum)
+ xp(ilocnum) = coord(1,iglob)
+ zp(ilocnum) = coord(2,iglob)
+ enddo
+
+ ! sorts (lexicographically?) ibool_interfaces and updates value
+ ! of total number of points nibool_interfaces_true(iinterface)
+ call sort_array_coordinates(num_nibool,xp,zp, &
+ ibool_dummy, &
+ reorder_interface,locval,ifseg, &
+ nibool_interfaces_true(iinterface), &
+ ind,ninseg,iwork,work)
+
+ ! checks that number of MPI points are still the same
+ num_points1 = num_points1 + num_nibool
+ num_points2 = num_points2 + nibool_interfaces_true(iinterface)
+ if( num_points1 /= num_points2 ) then
+ write(IOUT,*) 'error sorting MPI interface points:',myrank
+ write(IOUT,*) ' domain:',idomain
+ write(IOUT,*) ' interface:',iinterface,num_points1,num_points2
+ call exit_MPI('error sorting MPI interface')
+ endif
+
+ ! stores new order of ibool array
+ if( idomain == 1 ) then
+ ibool_interfaces_acoustic(1:num_nibool,iinterface) = ibool_dummy(:)
+ elseif( idomain == 2 ) then
+ ibool_interfaces_elastic(1:num_nibool,iinterface) = ibool_dummy(:)
+ elseif( idomain == 3 ) then
+ ibool_interfaces_poroelastic(1:num_nibool,iinterface) = ibool_dummy(:)
+ endif
+
+ ! cleanup temporary arrays
+ deallocate(xp)
+ deallocate(zp)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(reorder_interface)
+ deallocate(ibool_dummy)
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iwork)
+ deallocate(work)
+ enddo
+ enddo
+
+ ! cleanup
+ deallocate(nibool_interfaces_true)
+
+ ! outputs total number of MPI interface points
+ call MPI_ALLREDUCE(num_points2, num_points1, 1, MPI_INTEGER, &
+ MPI_SUM, MPI_COMM_WORLD, ier)
+ if( myrank == 0 .and. ipass == 1 ) then
+ write(IOUT,*) 'total MPI interface points: ',num_points1
+ endif
+
+ ! checks interfaces in acoustic domains
+ if ( ninterface_acoustic > 0) then
+
+ ! checks with assembly of test fields
+ allocate(test_flag_cr(npoin))
+ test_flag_cr(:) = 0._CUSTOM_REAL
+ count = 0
+ do ispec = 1, nspec
+ ! sets flags on global points
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ ! global index
+ iglob = ibool(i,j,ispec)
+
+ ! counts number of unique global points to set
+ if( nint(test_flag_cr(iglob)) == 0 ) count = count+1
+
+ ! sets identifier
+ test_flag_cr(iglob) = myrank + 1.0
+ enddo
+ enddo
+ enddo
+
+ max_nibool_interfaces = maxval(nibool_interfaces_acoustic(:))
+
+ allocate(tab_requests_send_recv_acoustic(ninterface_acoustic*2))
+ allocate(buffer_send_faces_vector_ac(max_nibool_interfaces,ninterface_acoustic))
+ allocate(buffer_recv_faces_vector_ac(max_nibool_interfaces,ninterface_acoustic))
+
+ inum = 0
+ do iinterface = 1, ninterface
+ inum = inum + nibool_interfaces_acoustic(iinterface)
+ enddo
+
+ call MPI_ALLREDUCE(inum, num_points2, 1, MPI_INTEGER, &
+ MPI_SUM, MPI_COMM_WORLD, ier)
+
+ if( myrank == 0 .and. ipass == 1 ) then
+ write(IOUT,*) ' acoustic interface points: ',num_points2
+ endif
+
+ ! adds contributions from different partitions to flag arrays
+ ! custom_real arrays
+ call assemble_MPI_vector_ac(test_flag_cr,npoin, &
+ ninterface, ninterface_acoustic,inum_interfaces_acoustic, &
+ max_interface_size, max_nibool_interfaces,&
+ ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+ tab_requests_send_recv_acoustic,buffer_send_faces_vector_ac, &
+ buffer_recv_faces_vector_ac, my_neighbours)
+
+ ! checks number of interface points
+ i = 0
+ do iglob=1,npoin
+ ! only counts flags with MPI contributions
+ if( test_flag_cr(iglob) > myrank+1.0_CUSTOM_REAL ) i = i + 1
+ enddo
+ call MPI_ALLREDUCE(inum, iglob, 1, MPI_INTEGER, &
+ MPI_SUM, MPI_COMM_WORLD, ier)
+
+ if( myrank == 0 .and. ipass == 1 ) then
+ write(IOUT,*) ' assembled acoustic MPI interface points:',iglob
+ endif
+ if( num_points2 /= iglob ) then
+ print*,'error assembly:',myrank
+ print*,' count = ',count
+ print*,' inum = ',inum
+ print*,' i = ',i
+ print*,' total: ',num_points2,' not equal to assembled ',iglob
+ call exit_MPI('error acoustic MPI assembly')
+ endif
+ deallocate(tab_requests_send_recv_acoustic)
+ deallocate(buffer_send_faces_vector_ac)
+ deallocate(buffer_recv_faces_vector_ac)
+
+ deallocate(test_flag_cr)
+
+ endif
+
+
+
+ end subroutine get_MPI
+
+#endif
Modified: seismo/2D/SPECFEM2D/trunk/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/gmat01.f90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/gmat01.f90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -196,8 +196,10 @@
C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
- call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare,H_biot,C_biot,M_biot,mu_fr,phi, &
- tortuosity,density(1),density(2),eta_f,val4,f0,freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
+ call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare, &
+ H_biot,C_biot,M_biot,mu_fr,phi, &
+ tortuosity,density(1),density(2),eta_f, &
+ val4,f0,freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
porosity_array(n) = val2
tortuosity_array(n) = val3
@@ -238,7 +240,7 @@
else
porosity_array(n) = 1.d0
endif
- elseif(indic == 2) then
+ elseif (indic == 2) then
density_array(1,n) = density(1)
! dummy poroelastcoef values, trick to avoid floating invalid
poroelastcoef(1,1,n) = lambda
@@ -255,7 +257,7 @@
Qp_array(n) = 15
Qs_array(n) = 15
porosity_array(n) = 0.d0
- else
+ elseif (indic == 3) then
density_array(1,n) = density(1)
density_array(2,n) = density(2)
poroelastcoef(1,1,n) = lambda_s
@@ -289,7 +291,7 @@
endif
elseif(indic == 2) then ! elastic (anisotropic)
write(IOUT,400) n,density(1),c11,c13,c15,c33,c35,c55
- else
+ elseif(indic == 3) then
! material is poroelastic (solid/fluid)
write(iout,500) n,sqrt(cpIsquare),sqrt(cpIIsquare),sqrt(cssquare)
write(iout,600) density(1),poisson_s,lambda_s,mu_s,kappa_s,young_s
Modified: seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -788,7 +788,8 @@
if ( (tab_size_interfaces(num_interface) < tab_size_interfaces(num_interface+1)) .and. &
(i == iproc .or. j == iproc) ) then
my_interfaces(num_interface) = 1
- my_nb_interfaces(num_interface) = tab_size_interfaces(num_interface+1) - tab_size_interfaces(num_interface)
+ my_nb_interfaces(num_interface) = tab_size_interfaces(num_interface+1) &
+ - tab_size_interfaces(num_interface)
endif
num_interface = num_interface + 1
enddo
@@ -798,55 +799,62 @@
else
do i = 0, nparts-1
- do j = i+1, nparts-1
- if ( my_interfaces(num_interface) == 1 ) then
- if ( i == iproc ) then
- write(IIN_database,*) j, my_nb_interfaces(num_interface)
- else
- write(IIN_database,*) i, my_nb_interfaces(num_interface)
- endif
+ do j = i+1, nparts-1
+ if ( my_interfaces(num_interface) == 1 ) then
+ if ( i == iproc ) then
+ write(IIN_database,*) j, my_nb_interfaces(num_interface)
+ else
+ write(IIN_database,*) i, my_nb_interfaces(num_interface)
+ endif
- do k = tab_size_interfaces(num_interface), tab_size_interfaces(num_interface+1)-1
- if ( i == iproc ) then
- local_elmnt = glob2loc_elmnts(tab_interfaces(k*5+0))+1
- else
- local_elmnt = glob2loc_elmnts(tab_interfaces(k*5+1))+1
- endif
+ do k = tab_size_interfaces(num_interface), tab_size_interfaces(num_interface+1)-1
+ if ( i == iproc ) then
+ local_elmnt = glob2loc_elmnts(tab_interfaces(k*5+0))+1
+ else
+ local_elmnt = glob2loc_elmnts(tab_interfaces(k*5+1))+1
+ endif
- if ( tab_interfaces(k*5+2) == 1 ) then
- do l = glob2loc_nodes_nparts(tab_interfaces(k*5+3)), &
+ if ( tab_interfaces(k*5+2) == 1 ) then
+ ! common node (single point)
+ do l = glob2loc_nodes_nparts(tab_interfaces(k*5+3)), &
glob2loc_nodes_nparts(tab_interfaces(k*5+3)+1)-1
- if ( glob2loc_nodes_parts(l) == iproc ) then
- local_nodes(1) = glob2loc_nodes(l)+1
- endif
- enddo
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(1) = glob2loc_nodes(l)+1
+ endif
+ enddo
- write(IIN_database,*) local_elmnt, tab_interfaces(k*5+2), local_nodes(1), -1
- else
- if ( tab_interfaces(k*5+2) == 2 ) then
- do l = glob2loc_nodes_nparts(tab_interfaces(k*5+3)), &
+ write(IIN_database,*) local_elmnt, tab_interfaces(k*5+2), &
+ local_nodes(1), -1
+ else
+ if ( tab_interfaces(k*5+2) == 2 ) then
+ ! common edge (two nodes)
+ ! first node
+ do l = glob2loc_nodes_nparts(tab_interfaces(k*5+3)), &
glob2loc_nodes_nparts(tab_interfaces(k*5+3)+1)-1
- if ( glob2loc_nodes_parts(l) == iproc ) then
- local_nodes(1) = glob2loc_nodes(l)+1
- endif
- enddo
- do l = glob2loc_nodes_nparts(tab_interfaces(k*5+4)), &
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(1) = glob2loc_nodes(l)+1
+ endif
+ enddo
+ ! second node
+ do l = glob2loc_nodes_nparts(tab_interfaces(k*5+4)), &
glob2loc_nodes_nparts(tab_interfaces(k*5+4)+1)-1
- if ( glob2loc_nodes_parts(l) == iproc ) then
- local_nodes(2) = glob2loc_nodes(l)+1
- endif
- enddo
- write(IIN_database,*) local_elmnt, tab_interfaces(k*5+2), local_nodes(1), local_nodes(2)
- else
- write(IIN_database,*) "erreur_write_interface_", tab_interfaces(k*5+2)
- endif
- endif
- enddo
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(2) = glob2loc_nodes(l)+1
+ endif
+ enddo
- endif
+ write(IIN_database,*) local_elmnt, tab_interfaces(k*5+2), &
+ local_nodes(1), local_nodes(2)
+ else
+ write(IIN_database,*) "erreur_write_interface_", tab_interfaces(k*5+2)
+ endif
+ endif
+ enddo
- num_interface = num_interface + 1
- enddo
+ endif
+
+ num_interface = num_interface + 1
+ enddo
enddo
endif
Added: seismo/2D/SPECFEM2D/trunk/prepare_assemble_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/prepare_assemble_MPI.F90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/prepare_assemble_MPI.F90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -0,0 +1,338 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 6.1
+! ------------------------------
+!
+! Copyright Universite de Pau, CNRS and INRIA, France,
+! and Princeton University / California Institute of Technology, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+!
+! This file contains subroutines related to assembling (of the mass matrix, potential_dot_dot and
+! accel_elastic, accels_poroelastic, accelw_poroelastic).
+! These subroutines are for the most part not used in the sequential version.
+!
+
+#ifdef USE_MPI
+
+!-----------------------------------------------
+! Determines the points that are on the interfaces with other partitions, to help
+! build the communication buffers, and determines which elements are considered 'inner'
+! (no points in common with other partitions) and 'outer' (at least one point in common
+! with neighbouring partitions).
+! We have both acoustic and (poro)elastic buffers, for coupling between acoustic and (poro)elastic elements
+! led us to have two sets of communications.
+!-----------------------------------------------
+ subroutine prepare_assemble_MPI(nspec,ibool,knods, ngnod,npoin, elastic, poroelastic, &
+ ninterface, max_interface_size, &
+ my_nelmnts_neighbours, my_interfaces, &
+ ibool_interfaces_acoustic, ibool_interfaces_elastic, &
+ ibool_interfaces_poroelastic, &
+ nibool_interfaces_acoustic, nibool_interfaces_elastic, &
+ nibool_interfaces_poroelastic, &
+ inum_interfaces_acoustic, inum_interfaces_elastic, &
+ inum_interfaces_poroelastic, &
+ ninterface_acoustic, ninterface_elastic, ninterface_poroelastic, &
+ mask_ispec_inner_outer )
+
+ implicit none
+
+ include 'constants.h'
+
+ integer, intent(in) :: nspec, npoin, ngnod
+ logical, dimension(nspec), intent(in) :: elastic, poroelastic
+ integer, dimension(ngnod,nspec), intent(in) :: knods
+ integer, dimension(NGLLX,NGLLZ,nspec), intent(in) :: ibool
+
+ integer :: ninterface
+ integer :: max_interface_size
+ integer, dimension(ninterface) :: my_nelmnts_neighbours
+ integer, dimension(4,max_interface_size,ninterface) :: my_interfaces
+ integer, dimension(NGLLX*max_interface_size,ninterface) :: &
+ ibool_interfaces_acoustic,ibool_interfaces_elastic,ibool_interfaces_poroelastic
+ integer, dimension(ninterface) :: &
+ nibool_interfaces_acoustic,nibool_interfaces_elastic,nibool_interfaces_poroelastic
+
+ integer, dimension(ninterface), intent(out) :: &
+ inum_interfaces_acoustic, inum_interfaces_elastic, inum_interfaces_poroelastic
+ integer, intent(out) :: ninterface_acoustic, ninterface_elastic, ninterface_poroelastic
+
+ logical, dimension(nspec), intent(inout) :: mask_ispec_inner_outer
+
+ ! local parameters
+ integer :: num_interface
+ integer :: ispec_interface
+ logical, dimension(npoin) :: mask_ibool_acoustic
+ logical, dimension(npoin) :: mask_ibool_elastic
+ logical, dimension(npoin) :: mask_ibool_poroelastic
+ integer :: ixmin, ixmax, izmin, izmax, ix, iz
+ integer, dimension(ngnod) :: n
+ integer :: e1, e2, itype, ispec, k, sens, iglob
+ integer :: npoin_interface_acoustic
+ integer :: npoin_interface_elastic
+ integer :: npoin_interface_poroelastic
+
+ ! initializes
+ ibool_interfaces_acoustic(:,:) = 0
+ nibool_interfaces_acoustic(:) = 0
+ ibool_interfaces_elastic(:,:) = 0
+ nibool_interfaces_elastic(:) = 0
+ ibool_interfaces_poroelastic(:,:) = 0
+ nibool_interfaces_poroelastic(:) = 0
+
+ do num_interface = 1, ninterface
+ npoin_interface_acoustic = 0
+ npoin_interface_elastic = 0
+ npoin_interface_poroelastic = 0
+ mask_ibool_acoustic(:) = .false.
+ mask_ibool_elastic(:) = .false.
+ mask_ibool_poroelastic(:) = .false.
+
+ do ispec_interface = 1, my_nelmnts_neighbours(num_interface)
+ ! element id
+ ispec = my_interfaces(1,ispec_interface,num_interface)
+ ! type of interface: 1 = common point, 2 = common edge
+ itype = my_interfaces(2,ispec_interface,num_interface)
+ ! element control node ids
+ do k = 1, ngnod
+ n(k) = knods(k,ispec)
+ end do
+ ! common node ids
+ e1 = my_interfaces(3,ispec_interface,num_interface)
+ e2 = my_interfaces(4,ispec_interface,num_interface)
+
+ call get_edge(ngnod, n, itype, e1, e2, ixmin, ixmax, izmin, izmax, sens)
+
+ do iz = izmin, izmax, sens
+ do ix = ixmin, ixmax, sens
+ ! global index
+ iglob = ibool(ix,iz,ispec)
+
+ ! checks to which material this common interface belongs
+ if ( elastic(ispec) ) then
+ ! elastic element
+ if(.not. mask_ibool_elastic(iglob)) then
+ mask_ibool_elastic(iglob) = .true.
+ npoin_interface_elastic = npoin_interface_elastic + 1
+ ibool_interfaces_elastic(npoin_interface_elastic,num_interface) = iglob
+ end if
+ else if ( poroelastic(ispec) ) then
+ ! poroelastic element
+ if(.not. mask_ibool_poroelastic(iglob)) then
+ mask_ibool_poroelastic(iglob) = .true.
+ npoin_interface_poroelastic = npoin_interface_poroelastic + 1
+ ibool_interfaces_poroelastic(npoin_interface_poroelastic,num_interface) = iglob
+ end if
+ else
+ ! acoustic element
+ if(.not. mask_ibool_acoustic(iglob)) then
+ mask_ibool_acoustic(iglob) = .true.
+ npoin_interface_acoustic = npoin_interface_acoustic + 1
+ ibool_interfaces_acoustic(npoin_interface_acoustic,num_interface) = iglob
+ end if
+ end if
+ end do
+ end do
+
+ end do
+
+ ! stores counters for interface points
+ nibool_interfaces_acoustic(num_interface) = npoin_interface_acoustic
+ nibool_interfaces_elastic(num_interface) = npoin_interface_elastic
+ nibool_interfaces_poroelastic(num_interface) = npoin_interface_poroelastic
+
+ ! sets inner/outer element flags
+ do ispec = 1, nspec
+ do iz = 1, NGLLZ
+ do ix = 1, NGLLX
+ if ( mask_ibool_acoustic(ibool(ix,iz,ispec)) &
+ .or. mask_ibool_elastic(ibool(ix,iz,ispec)) &
+ .or. mask_ibool_poroelastic(ibool(ix,iz,ispec)) ) then
+ mask_ispec_inner_outer(ispec) = .true.
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ end do
+
+ ! sets number of interfaces for each material domain
+ ninterface_acoustic = 0
+ ninterface_elastic = 0
+ ninterface_poroelastic = 0
+
+ do num_interface = 1, ninterface
+ ! acoustic
+ if ( nibool_interfaces_acoustic(num_interface) > 0 ) then
+ ninterface_acoustic = ninterface_acoustic + 1
+ inum_interfaces_acoustic(ninterface_acoustic) = num_interface
+ end if
+ ! elastic
+ if ( nibool_interfaces_elastic(num_interface) > 0 ) then
+ ninterface_elastic = ninterface_elastic + 1
+ inum_interfaces_elastic(ninterface_elastic) = num_interface
+ end if
+ ! poroelastic
+ if ( nibool_interfaces_poroelastic(num_interface) > 0 ) then
+ ninterface_poroelastic = ninterface_poroelastic + 1
+ inum_interfaces_poroelastic(ninterface_poroelastic) = num_interface
+ end if
+ end do
+
+ end subroutine prepare_assemble_MPI
+
+
+!-----------------------------------------------
+! Get the points (ixmin, ixmax, izmin and izmax) on an node/edge for one element.
+! 'sens' is used to have DO loops with increment equal to 'sens' (-/+1).
+!-----------------------------------------------
+ subroutine get_edge ( ngnod, n, itype, e1, e2, ixmin, ixmax, izmin, izmax, sens )
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: ngnod
+ integer, dimension(ngnod), intent(in) :: n
+ integer, intent(in) :: itype, e1, e2
+ integer, intent(out) :: ixmin, ixmax, izmin, izmax
+ integer, intent(out) :: sens
+
+ if ( itype == 1 ) then
+
+ ! common single point
+
+ ! checks which corner point is given
+ if ( e1 == n(1) ) then
+ ixmin = 1
+ ixmax = 1
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(2) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(3) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ if ( e1 == n(4) ) then
+ ixmin = 1
+ ixmax = 1
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ sens = 1
+
+ else if( itype == 2 ) then
+
+ ! common edge
+
+ ! checks which edge and corner points are given
+ if ( e1 == n(1) ) then
+ ixmin = 1
+ izmin = 1
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ izmax = 1
+ sens = 1
+ end if
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ izmax = NGLLZ
+ sens = 1
+ end if
+ end if
+ if ( e1 == n(2) ) then
+ ixmin = NGLLX
+ izmin = 1
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ izmax = NGLLZ
+ sens = 1
+ end if
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ izmax = 1
+ sens = -1
+ end if
+ end if
+ if ( e1 == n(3) ) then
+ ixmin = NGLLX
+ izmin = NGLLZ
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ izmax = NGLLZ
+ sens = -1
+ end if
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ izmax = 1
+ sens = -1
+ end if
+ end if
+ if ( e1 == n(4) ) then
+ ixmin = 1
+ izmin = NGLLZ
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ izmax = 1
+ sens = -1
+ end if
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ izmax = NGLLZ
+ sens = 1
+ end if
+ end if
+
+ else
+
+ call exit_MPI('ERROR get_edge unknown type')
+
+ end if
+
+ end subroutine get_edge
+
+#endif
Modified: seismo/2D/SPECFEM2D/trunk/prepare_color_image.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/prepare_color_image.F90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/prepare_color_image.F90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -273,14 +273,16 @@
double precision :: afactor,bfactor,cfactor,D_biot,H_biot,C_biot,&
M_biot,B_biot,cpIsquare,cpIIsquare,cssquare
double precision :: gamma1,gamma2,gamma3,gamma4,ratio
- integer :: i,j,ispec
+ integer :: i,j,k,ispec
#ifdef USE_MPI
double precision, dimension(:), allocatable :: data_pixel_recv
double precision, dimension(:), allocatable :: data_pixel_send
integer, dimension(:,:), allocatable :: num_pixel_recv
integer, dimension(:), allocatable :: nb_pixel_per_proc
integer, dimension(MPI_STATUS_SIZE) :: request_mpi_status
- integer :: ier,k,iproc
+ integer :: ier,iproc
+#else
+ integer :: dummy
#endif
! to display the P-velocity model in background on color images
@@ -425,8 +427,10 @@
deallocate(num_pixel_recv)
deallocate(data_pixel_recv)
endif
-
+#else
+ ! to avoid compiler warnings
+ dummy = myrank
+ dummy = nproc
#endif
-
-
+
end subroutine prepare_color_image_vp
Modified: seismo/2D/SPECFEM2D/trunk/read_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_databases.f90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/read_databases.f90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -413,9 +413,12 @@
allocate(knods_read(ngnod))
n = 0
do ispec = 1,nspec
+ ! format: #element_id #material_id #node_id1 #node_id2 #...
read(IIN,*) n,kmato_read,(knods_read(k), k=1,ngnod)
if(ipass == 1) then
+ ! material association
kmato(n) = kmato_read
+ ! element control node indices
knods(:,n)= knods_read(:)
else if(ipass == 2) then
kmato(perm(antecedent_list(n))) = kmato_read
@@ -474,14 +477,25 @@
integer :: num_interface,ie,my_interfaces_read
! initializes
- my_neighbours(:) = 0
+ my_neighbours(:) = -1
my_nelmnts_neighbours(:) = 0
- my_interfaces(:,:,:) = 0
+ my_interfaces(:,:,:) = -1
! reads in interfaces
do num_interface = 1, ninterface
+ ! format: #process_interface_id #number_of_elements_on_interface
+ ! where
+ ! process_interface_id = rank of (neighbor) process to share MPI interface with
+ ! number_of_elements_on_interface = number of interface elements
read(IIN,*) my_neighbours(num_interface), my_nelmnts_neighbours(num_interface)
+
+ ! loops over interface elements
do ie = 1, my_nelmnts_neighbours(num_interface)
+ ! format: #(1)spectral_element_id #(2)interface_type #(3)node_id1 #(4)node_id2
+ !
+ ! interface types:
+ ! 1 - corner point only
+ ! 2 - element edge
read(IIN,*) my_interfaces_read, my_interfaces(2,ie,num_interface), &
my_interfaces(3,ie,num_interface), my_interfaces(4,ie,num_interface)
@@ -506,7 +520,9 @@
subroutine read_databases_absorbing(myrank,ipass,nelemabs,nspec,anyabs, &
ibegin_bottom,iend_bottom,jbegin_right,jend_right, &
ibegin_top,iend_top,jbegin_left,jend_left, &
- numabs,codeabs,perm,antecedent_list)
+ numabs,codeabs,perm,antecedent_list, &
+ nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax, &
+ ib_right,ib_left,ib_bottom,ib_top)
! reads in absorbing edges
@@ -520,7 +536,10 @@
logical, dimension(4,nelemabs) :: codeabs
logical :: anyabs
integer, dimension(nspec) :: perm,antecedent_list
+ integer :: nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax
+ integer, dimension(nelemabs) :: ib_right,ib_left,ib_bottom,ib_top
+
! local parameters
integer :: inum,numabsread
logical :: codeabsread(4)
@@ -528,20 +547,33 @@
! initializes
codeabs(:,:) = .false.
+
ibegin_bottom(:) = 0
iend_bottom(:) = 0
ibegin_top(:) = 0
iend_top(:) = 0
+
jbegin_left(:) = 0
jend_left(:) = 0
jbegin_right(:) = 0
jend_right(:) = 0
+ nspec_xmin = 0
+ nspec_xmax = 0
+ nspec_zmin = 0
+ nspec_zmax = 0
+
+ ib_right(:) = 0
+ ib_left(:) = 0
+ ib_bottom(:) = 0
+ ib_top(:) = 0
+
! reads in absorbing edges
read(IIN,"(a80)") datlin
! reads in values
if( anyabs ) then
+ ! reads absorbing boundaries
do inum = 1,nelemabs
read(IIN,*) numabsread,codeabsread(1),codeabsread(2),codeabsread(3),&
codeabsread(4), ibegin_bottom(inum), iend_bottom(inum), &
@@ -565,9 +597,34 @@
codeabs(ILEFT,inum) = codeabsread(4)
enddo
+ ! boundary element numbering
+ do inum = 1,nelemabs
+ if (codeabs(IBOTTOM,inum)) then
+ nspec_zmin = nspec_zmin + 1
+ ib_bottom(inum) = nspec_zmin
+ endif
+ if (codeabs(IRIGHT,inum)) then
+ nspec_xmax = nspec_xmax + 1
+ ib_right(inum) = nspec_xmax
+ endif
+ if (codeabs(ITOP,inum)) then
+ nspec_zmax = nspec_zmax + 1
+ ib_top(inum) = nspec_zmax
+ endif
+ if (codeabs(ILEFT,inum)) then
+ nspec_xmin = nspec_xmin + 1
+ ib_left(inum) = nspec_xmin
+ endif
+ enddo
+
if (myrank == 0 .and. ipass == 1) then
write(IOUT,*)
write(IOUT,*) 'Number of absorbing elements: ',nelemabs
+ write(IOUT,*) ' nspec_xmin = ',nspec_xmin
+ write(IOUT,*) ' nspec_xmax = ',nspec_xmax
+ write(IOUT,*) ' nspec_zmin = ',nspec_zmin
+ write(IOUT,*) ' nspec_zmax = ',nspec_zmax
+ write(IOUT,*)
endif
endif
Modified: seismo/2D/SPECFEM2D/trunk/read_external_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_external_model.f90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/read_external_model.f90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -44,131 +44,145 @@
!========================================================================
-subroutine read_external_model(any_acoustic,any_elastic,any_poroelastic, &
+ subroutine read_external_model(any_acoustic,any_elastic,any_poroelastic, &
elastic,poroelastic,anisotropic,nspec,npoin,N_SLS,ibool, &
f0_attenuation,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent, &
inv_tau_sigma_nu1,inv_tau_sigma_nu2,phi_nu1,phi_nu2,Mu_nu1,Mu_nu2,&
coord,kmato,myrank,rhoext,vpext,vsext, &
- Qp_attenuationext,Qs_attenuationext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,READ_EXTERNAL_SEP_FILE)
+ Qp_attenuationext,Qs_attenuationext, &
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,READ_EXTERNAL_SEP_FILE)
implicit none
include "constants.h"
-integer :: nspec,myrank,npoin
-double precision :: f0_attenuation
+ integer :: nspec,myrank,npoin
+ double precision :: f0_attenuation
-! Mesh
-integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
-double precision, dimension(NDIM,npoin) :: coord
+ ! Mesh
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+ double precision, dimension(NDIM,npoin) :: coord
-! Material properties
-logical :: any_acoustic,any_elastic,any_poroelastic,READ_EXTERNAL_SEP_FILE
-integer, dimension(nspec) :: kmato
-logical, dimension(nspec) :: elastic,poroelastic
-double precision, dimension(NGLLX,NGLLZ,nspec) :: rhoext,vpext,vsext
+ ! Material properties
+ logical :: any_acoustic,any_elastic,any_poroelastic,READ_EXTERNAL_SEP_FILE
+ integer, dimension(nspec) :: kmato
+ logical, dimension(nspec) :: elastic,poroelastic
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: rhoext,vpext,vsext
-! for attenuation
-integer :: N_SLS
-double precision :: Mu_nu1_sent,Mu_nu2_sent
-double precision, dimension(N_SLS) :: inv_tau_sigma_nu1_sent,phi_nu1_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent
-double precision, dimension(NGLLX,NGLLZ,nspec,N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
-double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
-double precision, dimension(NGLLX,NGLLZ,nspec) :: Qp_attenuationext,Qs_attenuationext
+ ! for attenuation
+ integer :: N_SLS
+ double precision :: Mu_nu1_sent,Mu_nu2_sent
+ double precision, dimension(N_SLS) :: inv_tau_sigma_nu1_sent,phi_nu1_sent, &
+ inv_tau_sigma_nu2_sent,phi_nu2_sent
+ double precision, dimension(NGLLX,NGLLZ,nspec,N_SLS) :: inv_tau_sigma_nu1,phi_nu1, &
+ inv_tau_sigma_nu2,phi_nu2
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: Qp_attenuationext,Qs_attenuationext
-! for anisotropy
-logical, dimension(nspec) :: anisotropic
-double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c13ext,c15ext,c33ext,c35ext,c55ext
+ ! for anisotropy
+ logical, dimension(nspec) :: anisotropic
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c13ext,c15ext,c33ext,c35ext,c55ext
-! Local variables
-integer :: i,j,ispec,iglob
-double precision :: previous_vsext
-double precision :: tmp1, tmp2,tmp3
+ ! Local variables
+ integer :: i,j,ispec,iglob
+ double precision :: previous_vsext
+ double precision :: tmp1, tmp2,tmp3
-if(READ_EXTERNAL_SEP_FILE) then
- write(IOUT,*)
- write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
- write(IOUT,*) 'Assigning external velocity and density model (elastic (no attenuation) and/or acoustic)...'
- write(IOUT,*) 'Read outside SEG model...'
- write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+ if(READ_EXTERNAL_SEP_FILE) then
+ write(IOUT,*)
+ write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+ write(IOUT,*) 'Assigning external velocity and density model (elastic (no attenuation) and/or acoustic)...'
+ write(IOUT,*) 'Read outside SEG model...'
+ write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
- open(unit=1001,file='DATA/model_velocity.dat_input',status='unknown')
- do ispec = 1,nspec
- do j = 1,NGLLZ
- do i = 1,NGLLX
- iglob = ibool(i,j,ispec)
- read(1001,*) tmp1,tmp2,tmp3,rhoext(i,j,ispec),vpext(i,j,ispec),vsext(i,j,ispec)
-! vsext(i,j,ispec)=0.0
- ! Qp, Qs : dummy values. If attenuation needed than the "read" line and model_velocity.dat_input
- ! need to be modified to provide Qp & Qs values
- Qp_attenuationext(i,j,ispec) = 10.d0
- Qs_attenuationext(i,j,ispec) = 10.d0
- end do
- end do
+ open(unit=1001,file='DATA/model_velocity.dat_input',status='unknown')
+ do ispec = 1,nspec
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec)
+ read(1001,*) tmp1,tmp2,tmp3,rhoext(i,j,ispec),vpext(i,j,ispec),vsext(i,j,ispec)
+ ! vsext(i,j,ispec)=0.0
+ ! Qp, Qs : dummy values. If attenuation needed than the "read" line and model_velocity.dat_input
+ ! need to be modified to provide Qp & Qs values
+ Qp_attenuationext(i,j,ispec) = 10.d0
+ Qs_attenuationext(i,j,ispec) = 10.d0
end do
- close(1001)
+ end do
+ end do
+ close(1001)
- else
- do ispec = 1,nspec
- do j = 1,NGLLZ
- do i = 1,NGLLX
+ else
+ do ispec = 1,nspec
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
- iglob = ibool(i,j,ispec)
- call define_external_model(coord(1,iglob),coord(2,iglob),kmato(ispec),myrank,&
- rhoext(i,j,ispec),vpext(i,j,ispec),vsext(i,j,ispec), &
- Qp_attenuationext(i,j,ispec),Qs_attenuationext(i,j,ispec),&
- c11ext(i,j,ispec),c13ext(i,j,ispec),c15ext(i,j,ispec),c33ext(i,j,ispec),c35ext(i,j,ispec),c55ext(i,j,ispec))
- if((c11ext(i,j,ispec) /= 0) .or. (c13ext(i,j,ispec) /= 0) .or. (c15ext(i,j,ispec) /= 0) .or. &
- (c33ext(i,j,ispec) /= 0) .or. (c35ext(i,j,ispec) /= 0) .or. (c55ext(i,j,ispec) /= 0)) then
- ! vp, vs : dummy values, trick to avoid floating point errors
- vpext(i,j,ispec) = 20.d0
- vsext(i,j,ispec) = 10.d0
- end if
- end do
- end do
- end do
- end if
+ iglob = ibool(i,j,ispec)
+ call define_external_model(coord(1,iglob),coord(2,iglob),kmato(ispec),myrank,&
+ rhoext(i,j,ispec),vpext(i,j,ispec),vsext(i,j,ispec), &
+ Qp_attenuationext(i,j,ispec),Qs_attenuationext(i,j,ispec),&
+ c11ext(i,j,ispec),c13ext(i,j,ispec),c15ext(i,j,ispec), &
+ c33ext(i,j,ispec),c35ext(i,j,ispec),c55ext(i,j,ispec))
+
+ if((c11ext(i,j,ispec) /= 0) .or. (c13ext(i,j,ispec) /= 0) .or. (c15ext(i,j,ispec) /= 0) .or. &
+ (c33ext(i,j,ispec) /= 0) .or. (c35ext(i,j,ispec) /= 0) .or. (c55ext(i,j,ispec) /= 0)) then
+ ! vp, vs : dummy values, trick to avoid floating point errors
+ vpext(i,j,ispec) = 20.d0
+ vsext(i,j,ispec) = 10.d0
+ end if
+ end do
+ end do
+ end do
+ end if
- any_acoustic = .false.
- any_elastic = .false.
- any_poroelastic = .false.
- do ispec = 1,nspec
- previous_vsext = -1.d0
- do j = 1,NGLLZ
- do i = 1,NGLLX
- iglob = ibool(i,j,ispec)
- if(.not. (i == 1 .and. j == 1) .and. &
- ((vsext(i,j,ispec) >= TINYVAL .and. previous_vsext < TINYVAL) .or. &
- (vsext(i,j,ispec) < TINYVAL .and. previous_vsext >= TINYVAL))) &
- call exit_MPI('external velocity model cannot be both fluid and solid inside the same spectral element')
- if((c11ext(i,j,ispec) /= 0) .or. (c13ext(i,j,ispec) /= 0) .or. (c15ext(i,j,ispec) /= 0) .or. &
- (c33ext(i,j,ispec) /= 0) .or. (c35ext(i,j,ispec) /= 0) .or. (c55ext(i,j,ispec) /= 0)) then
- anisotropic(ispec) = .true.
- poroelastic(ispec) = .false.
- elastic(ispec) = .true.
- any_elastic = .true.
- Qp_attenuationext(i,j,ispec) = 10.d0
- Qs_attenuationext(i,j,ispec) = 10.d0
- elseif(vsext(i,j,ispec) < TINYVAL) then
- elastic(ispec) = .false.
- poroelastic(ispec) = .false.
- any_acoustic = .true.
- else
- poroelastic(ispec) = .false.
- elastic(ispec) = .true.
- any_elastic = .true.
- endif
- call attenuation_model(N_SLS,Qp_attenuationext(i,j,ispec),Qs_attenuationext(i,j,ispec), &
- f0_attenuation,inv_tau_sigma_nu1_sent,phi_nu1_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent)
- inv_tau_sigma_nu1(i,j,ispec,:) = inv_tau_sigma_nu1_sent(:)
- phi_nu1(i,j,ispec,:) = phi_nu1_sent(:)
- inv_tau_sigma_nu2(i,j,ispec,:) = inv_tau_sigma_nu2_sent(:)
- phi_nu2(i,j,ispec,:) = phi_nu2_sent(:)
- Mu_nu1(i,j,ispec) = Mu_nu1_sent
- Mu_nu2(i,j,ispec) = Mu_nu2_sent
- previous_vsext = vsext(i,j,ispec)
- enddo
- enddo
+ ! initializes
+ any_acoustic = .false.
+ any_elastic = .false.
+ any_poroelastic = .false.
+
+ anisotropic(:) = .false.
+ elastic(:) = .false.
+ poroelastic(:) = .false.
+
+ do ispec = 1,nspec
+ previous_vsext = -1.d0
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec)
+ if(.not. (i == 1 .and. j == 1) .and. &
+ ((vsext(i,j,ispec) >= TINYVAL .and. previous_vsext < TINYVAL) .or. &
+ (vsext(i,j,ispec) < TINYVAL .and. previous_vsext >= TINYVAL))) &
+ call exit_MPI('external velocity model cannot be both fluid and solid inside the same spectral element')
+
+ if((c11ext(i,j,ispec) /= 0) .or. (c13ext(i,j,ispec) /= 0) .or. (c15ext(i,j,ispec) /= 0) .or. &
+ (c33ext(i,j,ispec) /= 0) .or. (c35ext(i,j,ispec) /= 0) .or. (c55ext(i,j,ispec) /= 0)) then
+ anisotropic(ispec) = .true.
+ poroelastic(ispec) = .false.
+ elastic(ispec) = .true.
+ any_elastic = .true.
+ Qp_attenuationext(i,j,ispec) = 10.d0
+ Qs_attenuationext(i,j,ispec) = 10.d0
+ elseif(vsext(i,j,ispec) < TINYVAL) then
+ elastic(ispec) = .false.
+ poroelastic(ispec) = .false.
+ any_acoustic = .true.
+ else
+ poroelastic(ispec) = .false.
+ elastic(ispec) = .true.
+ any_elastic = .true.
+ endif
+
+ call attenuation_model(N_SLS,Qp_attenuationext(i,j,ispec),Qs_attenuationext(i,j,ispec), &
+ f0_attenuation,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
+ inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent)
+ inv_tau_sigma_nu1(i,j,ispec,:) = inv_tau_sigma_nu1_sent(:)
+ phi_nu1(i,j,ispec,:) = phi_nu1_sent(:)
+ inv_tau_sigma_nu2(i,j,ispec,:) = inv_tau_sigma_nu2_sent(:)
+ phi_nu2(i,j,ispec,:) = phi_nu2_sent(:)
+ Mu_nu1(i,j,ispec) = Mu_nu1_sent
+ Mu_nu2(i,j,ispec) = Mu_nu2_sent
+ previous_vsext = vsext(i,j,ispec)
enddo
+ enddo
+ enddo
-end subroutine read_external_model
+ end subroutine read_external_model
Added: seismo/2D/SPECFEM2D/trunk/sort_array_coordinates.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/sort_array_coordinates.F90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/sort_array_coordinates.F90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -0,0 +1,241 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 6.1
+! ------------------------------
+!
+! Copyright Universite de Pau, CNRS and INRIA, France,
+! and Princeton University / California Institute of Technology, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+
+#ifdef USE_MPI
+
+! subroutines to sort MPI buffers to assemble between chunks
+
+ subroutine sort_array_coordinates(npointot,x,z,ibool,iglob,loc,ifseg, &
+ nglob,ind,ninseg,iwork,work)
+
+! this routine MUST be in double precision to avoid sensitivity
+! to roundoff errors in the coordinates of the points
+!
+! returns: sorted indexing array (ibool), reordering array (iglob) & number of global points (nglob)
+
+ implicit none
+
+ include "constants.h"
+
+ integer,intent(in) :: npointot
+ integer,intent(out) :: nglob
+
+ integer,intent(inout) :: ibool(npointot)
+
+ integer iglob(npointot),loc(npointot)
+ integer ind(npointot),ninseg(npointot)
+ logical ifseg(npointot)
+ double precision,intent(in) :: x(npointot),z(npointot)
+ integer iwork(npointot)
+ double precision work(npointot)
+
+ ! local parameters
+ integer ipoin,i,j
+ integer nseg,ioff,iseg,ig
+ ! define a tolerance, normalized radius is 1., so let's use a small value
+ double precision,parameter :: xtol = SMALLVALTOL
+
+ ! establish initial pointers
+ do ipoin=1,npointot
+ loc(ipoin)=ipoin
+ enddo
+
+ ifseg(:)=.false.
+
+ nseg=1
+ ifseg(1)=.true.
+ ninseg(1)=npointot
+
+ do j=1,NDIM
+
+ ! sort within each segment
+ ioff=1
+ do iseg=1,nseg
+
+ if(j == 1) then
+ call rank_buffers(x(ioff),ind,ninseg(iseg))
+ else if(j == 2) then
+ call rank_buffers(z(ioff),ind,ninseg(iseg))
+ endif
+
+ call swap_all_buffers(ibool(ioff),loc(ioff), &
+ x(ioff),z(ioff),iwork,work,ind,ninseg(iseg))
+
+ ioff=ioff+ninseg(iseg)
+ enddo
+
+ ! check for jumps in current coordinate
+ if(j == 1) then
+ do i=2,npointot
+ if(dabs(x(i)-x(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else if(j == 2) then
+ do i=2,npointot
+ if(dabs(z(i)-z(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ endif
+
+ ! count up number of different segments
+ nseg=0
+ do i=1,npointot
+ if(ifseg(i)) then
+ nseg=nseg+1
+ ninseg(nseg)=1
+ else
+ ninseg(nseg)=ninseg(nseg)+1
+ endif
+ enddo
+ enddo
+
+ ! assign global node numbers (now sorted lexicographically)
+ ig=0
+ do i=1,npointot
+ if(ifseg(i)) ig=ig+1
+ iglob(loc(i))=ig
+ enddo
+
+ nglob=ig
+
+ end subroutine sort_array_coordinates
+
+! -------------------- library for sorting routine ------------------
+
+! sorting routines put here in same file to allow for inlining
+
+ subroutine rank_buffers(A,IND,N)
+!
+! Use Heap Sort (Numerical Recipes)
+!
+ implicit none
+
+ integer n
+ double precision A(n)
+ integer IND(n)
+
+ integer i,j,l,ir,indx
+ double precision q
+
+ do j=1,n
+ IND(j)=j
+ enddo
+
+ if(n == 1) return
+
+ L=n/2+1
+ ir=n
+ 100 CONTINUE
+ IF(l>1) THEN
+ l=l-1
+ indx=ind(l)
+ q=a(indx)
+ ELSE
+ indx=ind(ir)
+ q=a(indx)
+ ind(ir)=ind(1)
+ ir=ir-1
+ if (ir == 1) then
+ ind(1)=indx
+ return
+ endif
+ ENDIF
+ i=l
+ j=l+l
+ 200 CONTINUE
+ IF(J <= IR) THEN
+ IF(J < IR) THEN
+ IF(A(IND(j)) < A(IND(j+1))) j=j+1
+ ENDIF
+ IF (q < A(IND(j))) THEN
+ IND(I)=IND(J)
+ I=J
+ J=J+J
+ ELSE
+ J=IR+1
+ ENDIF
+ goto 200
+ ENDIF
+ IND(I)=INDX
+ goto 100
+ end subroutine rank_buffers
+
+! -------------------------------------------------------------------
+
+ subroutine swap_all_buffers(IA,IB,A,B,IW,W,ind,n)
+!
+! swap arrays IA, IB, A and B according to addressing in array IND
+!
+ implicit none
+
+ integer n
+
+ integer IND(n)
+ integer IA(n),IB(n),IW(n)
+ double precision A(n),B(n),W(n)
+
+ integer i
+
+ do i=1,n
+ W(i)=A(i)
+ IW(i)=IA(i)
+ enddo
+
+ do i=1,n
+ A(i)=W(ind(i))
+ IA(i)=IW(ind(i))
+ enddo
+
+ do i=1,n
+ W(i)=B(i)
+ IW(i)=IB(i)
+ enddo
+
+ do i=1,n
+ B(i)=W(ind(i))
+ IB(i)=IW(ind(i))
+ enddo
+
+ end subroutine swap_all_buffers
+
+#endif
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2011-02-21 20:39:45 UTC (rev 17930)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2011-02-22 01:39:25 UTC (rev 17931)
@@ -919,25 +919,33 @@
!-------------------------------------------------------------------------------
!---- determine if each spectral element is elastic, poroelastic, or acoustic
!-------------------------------------------------------------------------------
+ ! initializes
any_acoustic = .false.
any_elastic = .false.
any_poroelastic = .false.
- all_anisotropic = .false.
+
anisotropic(:) = .false.
+ elastic(:) = .false.
+ poroelastic(:) = .false.
+
+ ! loops over all elements
do ispec = 1,nspec
- if(porosity(kmato(ispec)) == 1.d0) then ! acoustic domain
+ if( nint(porosity(kmato(ispec))) == 1 ) then
+ ! acoustic domain
elastic(ispec) = .false.
poroelastic(ispec) = .false.
any_acoustic = .true.
- elseif(porosity(kmato(ispec)) < TINYVAL) then ! elastic domain
+ elseif( porosity(kmato(ispec)) < TINYVAL) then
+ ! elastic domain
elastic(ispec) = .true.
poroelastic(ispec) = .false.
any_elastic = .true.
if(any(anisotropy(:,kmato(ispec)) /= 0)) then
anisotropic(ispec) = .true.
end if
- else ! poroelastic domain
+ else
+ ! poroelastic domain
elastic(ispec) = .false.
poroelastic(ispec) = .true.
any_poroelastic = .true.
@@ -1093,38 +1101,42 @@
call read_databases_absorbing(myrank,ipass,nelemabs,nspec,anyabs, &
ibegin_bottom,iend_bottom,jbegin_right,jend_right, &
ibegin_top,iend_top,jbegin_left,jend_left, &
- numabs,codeabs,perm,antecedent_list)
+ numabs,codeabs,perm,antecedent_list, &
+ nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax, &
+ ib_right,ib_left,ib_bottom,ib_top)
+
if( anyabs ) then
- nspec_xmin = 0
- nspec_xmax = 0
- nspec_zmin = 0
- nspec_zmax = 0
+! nspec_xmin = 0
+! nspec_xmax = 0
+! nspec_zmin = 0
+! nspec_zmax = 0
! if(ipass == 1) then
! allocate(ib_left(nelemabs))
! allocate(ib_right(nelemabs))
! allocate(ib_bottom(nelemabs))
! allocate(ib_top(nelemabs))
! endif
- do inum = 1,nelemabs
- if (codeabs(IBOTTOM,inum)) then
- nspec_zmin = nspec_zmin + 1
- ib_bottom(inum) = nspec_zmin
- endif
- if (codeabs(IRIGHT,inum)) then
- nspec_xmax = nspec_xmax + 1
- ib_right(inum) = nspec_xmax
- endif
- if (codeabs(ITOP,inum)) then
- nspec_zmax = nspec_zmax + 1
- ib_top(inum) = nspec_zmax
- endif
- if (codeabs(ILEFT,inum)) then
- nspec_xmin = nspec_xmin + 1
- ib_left(inum) = nspec_xmin
- endif
- enddo
+
+! do inum = 1,nelemabs
+! if (codeabs(IBOTTOM,inum)) then
+! nspec_zmin = nspec_zmin + 1
+! ib_bottom(inum) = nspec_zmin
+! endif
+! if (codeabs(IRIGHT,inum)) then
+! nspec_xmax = nspec_xmax + 1
+! ib_right(inum) = nspec_xmax
+! endif
+! if (codeabs(ITOP,inum)) then
+! nspec_zmax = nspec_zmax + 1
+! ib_top(inum) = nspec_zmax
+! endif
+! if (codeabs(ILEFT,inum)) then
+! nspec_xmin = nspec_xmin + 1
+! ib_left(inum) = nspec_xmin
+! endif
+! enddo
! Files to save absorbed waves needed to reconstruct backward wavefield for adjoint method
if(ipass == 1) then
@@ -1171,13 +1183,13 @@
endif
endif
- if (myrank == 0 ) then
- write(IOUT,*)
- write(IOUT,*) 'nspec_xmin = ',nspec_xmin
- write(IOUT,*) 'nspec_xmax = ',nspec_xmax
- write(IOUT,*) 'nspec_zmin = ',nspec_zmin
- write(IOUT,*) 'nspec_zmax = ',nspec_zmax
- endif
+! if (myrank == 0 ) then
+! write(IOUT,*)
+! write(IOUT,*) 'nspec_xmin = ',nspec_xmin
+! write(IOUT,*) 'nspec_xmax = ',nspec_xmax
+! write(IOUT,*) 'nspec_zmin = ',nspec_zmin
+! write(IOUT,*) 'nspec_zmax = ',nspec_zmax
+! endif
else
@@ -1803,32 +1815,39 @@
c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,READ_EXTERNAL_SEP_FILE)
end if
+!
+!---- perform basic checks on parameters read
+!
+ all_anisotropic = .false.
if(count(anisotropic(:) .eqv. .true.) == nspec) all_anisotropic = .true.
+
if(all_anisotropic .and. anyabs) &
- stop 'Cannot put absorbing boundaries if anisotropic materials along edges'
+ call exit_MPI('Cannot put absorbing boundaries if anisotropic materials along edges')
+
if(TURN_ATTENUATION_ON .and. all_anisotropic) then
- stop 'Cannot turn attenuation on in anisotropic materials'
+ call exit_MPI('Cannot turn attenuation on in anisotropic materials')
end if
-!
-!---- perform basic checks on parameters read
-!
+ ! global domain flags
any_elastic_glob = any_elastic
#ifdef USE_MPI
- call MPI_ALLREDUCE(any_elastic, any_elastic_glob, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE(any_elastic, any_elastic_glob, 1, MPI_LOGICAL, &
+ MPI_LOR, MPI_COMM_WORLD, ier)
#endif
any_poroelastic_glob = any_poroelastic
#ifdef USE_MPI
- call MPI_ALLREDUCE(any_poroelastic, any_poroelastic_glob, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE(any_poroelastic, any_poroelastic_glob, 1, MPI_LOGICAL, &
+ MPI_LOR, MPI_COMM_WORLD, ier)
#endif
any_acoustic_glob = any_acoustic
#ifdef USE_MPI
- call MPI_ALLREDUCE(any_acoustic, any_acoustic_glob, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE(any_acoustic, any_acoustic_glob, 1, MPI_LOGICAL, &
+ MPI_LOR, MPI_COMM_WORLD, ier)
#endif
-! for acoustic
+ ! for acoustic
if(TURN_ATTENUATION_ON .and. .not. any_elastic_glob) &
call exit_MPI('currently cannot have attenuation if acoustic/poroelastic simulation only')
@@ -2420,17 +2439,25 @@
#ifdef USE_MPI
if ( nproc > 1 ) then
-! preparing for MPI communications
+
+ ! preparing for MPI communications
if(ipass == 1) allocate(mask_ispec_inner_outer(nspec))
mask_ispec_inner_outer(:) = .false.
- call prepare_assemble_MPI (nspec,ibool,knods,ngnod,npoin,elastic,poroelastic, &
- ninterface, max_interface_size,my_nelmnts_neighbours, my_interfaces, &
- ibool_interfaces_acoustic, ibool_interfaces_elastic, ibool_interfaces_poroelastic, &
- nibool_interfaces_acoustic, nibool_interfaces_elastic, nibool_interfaces_poroelastic, &
- inum_interfaces_acoustic, inum_interfaces_elastic, inum_interfaces_poroelastic, &
- ninterface_acoustic, ninterface_elastic, ninterface_poroelastic,mask_ispec_inner_outer)
+ call get_MPI(nspec,ibool,knods,ngnod,npoin,elastic,poroelastic, &
+ ninterface, max_interface_size, &
+ my_nelmnts_neighbours,my_interfaces,my_neighbours, &
+ ibool_interfaces_acoustic, ibool_interfaces_elastic, &
+ ibool_interfaces_poroelastic, &
+ nibool_interfaces_acoustic, nibool_interfaces_elastic, &
+ nibool_interfaces_poroelastic, &
+ inum_interfaces_acoustic, inum_interfaces_elastic, &
+ inum_interfaces_poroelastic, &
+ ninterface_acoustic, ninterface_elastic, ninterface_poroelastic, &
+ mask_ispec_inner_outer, &
+ myrank,ipass,coord)
+
nspec_outer = count(mask_ispec_inner_outer)
nspec_inner = nspec - nspec_outer
@@ -2439,7 +2466,7 @@
allocate(ispec_inner_to_glob(nspec_inner))
endif
-! building of corresponding arrays between inner/outer elements and their global number
+ ! building of corresponding arrays between inner/outer elements and their global number
if(ipass == 1) then
num_ispec_outer = 0
num_ispec_inner = 0
@@ -2454,6 +2481,7 @@
enddo
endif
+ ! buffers for MPI communications
max_ibool_interfaces_size_ac = maxval(nibool_interfaces_acoustic(:))
max_ibool_interfaces_size_el = 3*maxval(nibool_interfaces_elastic(:))
max_ibool_interfaces_size_po = NDIM*maxval(nibool_interfaces_poroelastic(:))
@@ -2472,12 +2500,17 @@
endif
! assembling the mass matrix
- call assemble_MPI_scalar(rmass_inverse_acoustic,rmass_inverse_elastic,rmass_s_inverse_poroelastic, &
- rmass_w_inverse_poroelastic,npoin, &
- ninterface, max_interface_size, max_ibool_interfaces_size_ac, max_ibool_interfaces_size_el, &
- max_ibool_interfaces_size_po, &
- ibool_interfaces_acoustic,ibool_interfaces_elastic,ibool_interfaces_poroelastic, &
- nibool_interfaces_acoustic,nibool_interfaces_elastic,nibool_interfaces_poroelastic,my_neighbours)
+ call assemble_MPI_scalar(rmass_inverse_acoustic, &
+ rmass_inverse_elastic, &
+ rmass_s_inverse_poroelastic, &
+ rmass_w_inverse_poroelastic,npoin, &
+ ninterface, max_interface_size, max_ibool_interfaces_size_ac, &
+ max_ibool_interfaces_size_el, &
+ max_ibool_interfaces_size_po, &
+ ibool_interfaces_acoustic,ibool_interfaces_elastic, &
+ ibool_interfaces_poroelastic, &
+ nibool_interfaces_acoustic,nibool_interfaces_elastic, &
+ nibool_interfaces_poroelastic,my_neighbours)
else
ninterface_acoustic = 0
@@ -3893,104 +3926,119 @@
! update displacement using finite-difference time scheme (Newmark)
if(any_elastic) then
- displ_elastic = displ_elastic + deltat*veloc_elastic + deltatsquareover2*accel_elastic
+ displ_elastic = displ_elastic &
+ + deltat*veloc_elastic &
+ + deltatsquareover2*accel_elastic
veloc_elastic = veloc_elastic + deltatover2*accel_elastic
accel_elastic = ZERO
- if(SIMULATION_TYPE == 2) then ! Adjoint calculation
- b_displ_elastic = b_displ_elastic + b_deltat*b_veloc_elastic + b_deltatsquareover2*b_accel_elastic
- b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic
- b_accel_elastic = ZERO
- endif
+ if(SIMULATION_TYPE == 2) then ! Adjoint calculation
+ b_displ_elastic = b_displ_elastic &
+ + b_deltat*b_veloc_elastic &
+ + b_deltatsquareover2*b_accel_elastic
+ b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic
+ b_accel_elastic = ZERO
+ endif
endif
if(any_poroelastic) then
-!for the solid
- displs_poroelastic = displs_poroelastic + deltat*velocs_poroelastic + deltatsquareover2*accels_poroelastic
+ !for the solid
+ displs_poroelastic = displs_poroelastic &
+ + deltat*velocs_poroelastic &
+ + deltatsquareover2*accels_poroelastic
velocs_poroelastic = velocs_poroelastic + deltatover2*accels_poroelastic
accels_poroelastic = ZERO
-!for the fluid
- displw_poroelastic = displw_poroelastic + deltat*velocw_poroelastic + deltatsquareover2*accelw_poroelastic
+ !for the fluid
+ displw_poroelastic = displw_poroelastic &
+ + deltat*velocw_poroelastic &
+ + deltatsquareover2*accelw_poroelastic
velocw_poroelastic = velocw_poroelastic + deltatover2*accelw_poroelastic
accelw_poroelastic = ZERO
- if(SIMULATION_TYPE == 2) then ! Adjoint calculation
-!for the solid
- b_displs_poroelastic = b_displs_poroelastic + b_deltat*b_velocs_poroelastic + b_deltatsquareover2*b_accels_poroelastic
- b_velocs_poroelastic = b_velocs_poroelastic + b_deltatover2*b_accels_poroelastic
- b_accels_poroelastic = ZERO
-!for the fluid
- b_displw_poroelastic = b_displw_poroelastic + b_deltat*b_velocw_poroelastic + b_deltatsquareover2*b_accelw_poroelastic
- b_velocw_poroelastic = b_velocw_poroelastic + b_deltatover2*b_accelw_poroelastic
- b_accelw_poroelastic = ZERO
- endif
+ if(SIMULATION_TYPE == 2) then ! Adjoint calculation
+ !for the solid
+ b_displs_poroelastic = b_displs_poroelastic &
+ + b_deltat*b_velocs_poroelastic &
+ + b_deltatsquareover2*b_accels_poroelastic
+ b_velocs_poroelastic = b_velocs_poroelastic + b_deltatover2*b_accels_poroelastic
+ b_accels_poroelastic = ZERO
+ !for the fluid
+ b_displw_poroelastic = b_displw_poroelastic &
+ + b_deltat*b_velocw_poroelastic &
+ + b_deltatsquareover2*b_accelw_poroelastic
+ b_velocw_poroelastic = b_velocw_poroelastic + b_deltatover2*b_accelw_poroelastic
+ b_accelw_poroelastic = ZERO
+ endif
endif
!--------------------------------------------------------------------------------------------
! implement viscous attenuation for poroelastic media
!
- if(TURN_VISCATTENUATION_ON .and. any_poroelastic) then
+ if(TURN_VISCATTENUATION_ON .and. any_poroelastic) then
! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
! loop over spectral elements
- do ispec = 1,nspec
+ do ispec = 1,nspec
- etal_f = poroelastcoef(2,2,kmato(ispec))
- permlxx = permeability(1,kmato(ispec))
- permlxz = permeability(2,kmato(ispec))
- permlzz = permeability(3,kmato(ispec))
+ etal_f = poroelastcoef(2,2,kmato(ispec))
+ permlxx = permeability(1,kmato(ispec))
+ permlxz = permeability(2,kmato(ispec))
+ permlzz = permeability(3,kmato(ispec))
- ! calcul of the inverse of k
+ ! calcul of the inverse of k
- detk = permlxx*permlzz - permlxz*permlxz
+ detk = permlxx*permlzz - permlxz*permlxz
- if(detk /= ZERO) then
- invpermlxx = permlzz/detk
- invpermlxz = -permlxz/detk
- invpermlzz = permlxx/detk
- else
- stop 'Permeability matrix is not invertible'
- endif
+ if(detk /= ZERO) then
+ invpermlxx = permlzz/detk
+ invpermlxz = -permlxz/detk
+ invpermlzz = permlxx/detk
+ else
+ stop 'Permeability matrix is not invertible'
+ endif
-! relaxed viscous coef
- bl_relaxed(1) = etal_f*invpermlxx
- bl_relaxed(2) = etal_f*invpermlxz
- bl_relaxed(3) = etal_f*invpermlzz
+ ! relaxed viscous coef
+ bl_relaxed(1) = etal_f*invpermlxx
+ bl_relaxed(2) = etal_f*invpermlxz
+ bl_relaxed(3) = etal_f*invpermlzz
- do j=1,NGLLZ
- do i=1,NGLLX
+ do j=1,NGLLZ
+ do i=1,NGLLX
- iglob = ibool(i,j,ispec)
+ iglob = ibool(i,j,ispec)
- viscox_loc(i,j) = velocw_poroelastic(1,iglob)*bl_relaxed(1) + &
+ viscox_loc(i,j) = velocw_poroelastic(1,iglob)*bl_relaxed(1) + &
velocw_poroelastic(2,iglob)*bl_relaxed(2)
- viscoz_loc(i,j) = velocw_poroelastic(1,iglob)*bl_relaxed(2) + &
+ viscoz_loc(i,j) = velocw_poroelastic(1,iglob)*bl_relaxed(2) + &
velocw_poroelastic(2,iglob)*bl_relaxed(3)
-! evolution rx_viscous
- Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscox(i,j,ispec)
- Snp1 = - (1.d0 - theta_e/theta_s)/theta_s*viscox_loc(i,j)
- rx_viscous(i,j,ispec) = alphaval * rx_viscous(i,j,ispec) + betaval * Sn + gammaval * Snp1
+ ! evolution rx_viscous
+ Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscox(i,j,ispec)
+ Snp1 = - (1.d0 - theta_e/theta_s)/theta_s*viscox_loc(i,j)
+ rx_viscous(i,j,ispec) = alphaval * rx_viscous(i,j,ispec) &
+ + betaval * Sn + gammaval * Snp1
-! evolution rz_viscous
- Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscoz(i,j,ispec)
- Snp1 = - (1.d0 - theta_e/theta_s)/theta_s*viscoz_loc(i,j)
- rz_viscous(i,j,ispec) = alphaval * rz_viscous(i,j,ispec) + betaval * Sn + gammaval * Snp1
+ ! evolution rz_viscous
+ Sn = - (1.d0 - theta_e/theta_s)/theta_s*viscoz(i,j,ispec)
+ Snp1 = - (1.d0 - theta_e/theta_s)/theta_s*viscoz_loc(i,j)
+ rz_viscous(i,j,ispec) = alphaval * rz_viscous(i,j,ispec) &
+ + betaval * Sn + gammaval * Snp1
- enddo
- enddo
+ enddo
+ enddo
-! save visco for Runge-Kutta scheme
- viscox(:,:,ispec) = viscox_loc(:,:)
- viscoz(:,:,ispec) = viscoz_loc(:,:)
+ ! save visco for Runge-Kutta scheme
+ viscox(:,:,ispec) = viscox_loc(:,:)
+ viscoz(:,:,ispec) = viscoz_loc(:,:)
- enddo ! end of spectral element loop
- endif ! end of viscous attenuation for porous media
+ enddo ! end of spectral element loop
+ endif ! end of viscous attenuation for porous media
!-----------------------------------------
if(any_acoustic) then
+ ! Newmark time scheme
potential_acoustic = potential_acoustic &
+ deltat*potential_dot_acoustic &
+ deltatsquareover2*potential_dot_dot_acoustic
@@ -4007,7 +4055,7 @@
b_potential_dot_dot_acoustic = ZERO
endif
-! free surface for an acoustic medium
+ ! free surface for an acoustic medium
if ( nelem_acoustic_surface > 0 ) then
call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,acoustic_surface, &
@@ -4069,10 +4117,9 @@
endif
-
+ ! stores absorbing boundary contributions into files
if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
-
-!--- left absorbing boundary
+ !--- left absorbing boundary
if(nspec_xmin >0) then
do ispec = 1,nspec_xmin
do i=1,NGLLZ
@@ -4080,43 +4127,39 @@
enddo
enddo
endif
+ !--- right absorbing boundary
+ if(nspec_xmax >0) then
+ do ispec = 1,nspec_xmax
+ do i=1,NGLLZ
+ write(66) b_absorb_acoustic_right(i,ispec,it)
+ enddo
+ enddo
+ endif
+ !--- bottom absorbing boundary
+ if(nspec_zmin >0) then
+ do ispec = 1,nspec_zmin
+ do i=1,NGLLX
+ write(67) b_absorb_acoustic_bottom(i,ispec,it)
+ enddo
+ enddo
+ endif
+ !--- top absorbing boundary
+ if(nspec_zmax >0) then
+ do ispec = 1,nspec_zmax
+ do i=1,NGLLX
+ write(68) b_absorb_acoustic_top(i,ispec,it)
+ enddo
+ enddo
+ endif
+ endif ! if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)
-!--- right absorbing boundary
- if(nspec_xmax >0) then
- do ispec = 1,nspec_xmax
- do i=1,NGLLZ
- write(66) b_absorb_acoustic_right(i,ispec,it)
- enddo
- enddo
- endif
+ endif ! end of test if any acoustic element
-!--- bottom absorbing boundary
- if(nspec_zmin >0) then
- do ispec = 1,nspec_zmin
- do i=1,NGLLX
- write(67) b_absorb_acoustic_bottom(i,ispec,it)
- enddo
- enddo
- endif
-
-!--- top absorbing boundary
- if(nspec_zmax >0) then
- do ispec = 1,nspec_zmax
- do i=1,NGLLX
- write(68) b_absorb_acoustic_top(i,ispec,it)
- enddo
- enddo
- endif
-
- endif ! if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)
-
- endif ! end of test if any acoustic element
-
! *********************************************************
! ************* add coupling with the elastic side
! *********************************************************
- if(coupled_acoustic_elastic) then
+ if(coupled_acoustic_elastic) then
! loop on all the coupling edges
do inum = 1,num_fluid_solid_edges
@@ -4141,8 +4184,8 @@
displ_z = displ_elastic(3,iglob)
if(SIMULATION_TYPE == 2) then
- b_displ_x = b_displ_elastic(1,iglob)
- b_displ_z = b_displ_elastic(3,iglob)
+ b_displ_x = b_displ_elastic(1,iglob)
+ b_displ_z = b_displ_elastic(3,iglob)
endif
! get point values for the acoustic side
@@ -4161,28 +4204,28 @@
jacobian1D = sqrt(xxi**2 + zxi**2)
nx = - zxi / jacobian1D
nz = + xxi / jacobian1D
- weight = jacobian1D * wxgll(i)
+ weight = jacobian1D * wxgll(i)
elseif(iedge_acoustic == IBOTTOM)then
xxi = + gammaz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
zxi = - gammax(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
jacobian1D = sqrt(xxi**2 + zxi**2)
nx = + zxi / jacobian1D
nz = - xxi / jacobian1D
- weight = jacobian1D * wxgll(i)
+ weight = jacobian1D * wxgll(i)
elseif(iedge_acoustic ==ILEFT)then
xgamma = - xiz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
zgamma = + xix(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
jacobian1D = sqrt(xgamma**2 + zgamma**2)
nx = - zgamma / jacobian1D
nz = + xgamma / jacobian1D
- weight = jacobian1D * wzgll(j)
+ weight = jacobian1D * wzgll(j)
elseif(iedge_acoustic ==IRIGHT)then
xgamma = - xiz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
zgamma = + xix(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
jacobian1D = sqrt(xgamma**2 + zgamma**2)
nx = + zgamma / jacobian1D
nz = - xgamma / jacobian1D
- weight = jacobian1D * wzgll(j)
+ weight = jacobian1D * wzgll(j)
endif
! compute dot product
@@ -4199,7 +4242,7 @@
enddo
- endif
+ endif
! *********************************************************
! ************* add coupling with the poroelastic side
@@ -4234,11 +4277,11 @@
displw_z = displw_poroelastic(2,iglob)
if(SIMULATION_TYPE == 2) then
- b_displ_x = b_displs_poroelastic(1,iglob)
- b_displ_z = b_displs_poroelastic(2,iglob)
+ b_displ_x = b_displs_poroelastic(1,iglob)
+ b_displ_z = b_displs_poroelastic(2,iglob)
- b_displw_x = b_displw_poroelastic(1,iglob)
- b_displw_z = b_displw_poroelastic(2,iglob)
+ b_displw_x = b_displw_poroelastic(1,iglob)
+ b_displw_z = b_displw_poroelastic(2,iglob)
endif
! get point values for the acoustic side
@@ -4258,28 +4301,28 @@
jacobian1D = sqrt(xxi**2 + zxi**2)
nx = - zxi / jacobian1D
nz = + xxi / jacobian1D
- weight = jacobian1D * wxgll(i)
+ weight = jacobian1D * wxgll(i)
elseif(iedge_acoustic == IBOTTOM)then
xxi = + gammaz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
zxi = - gammax(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
jacobian1D = sqrt(xxi**2 + zxi**2)
nx = + zxi / jacobian1D
nz = - xxi / jacobian1D
- weight = jacobian1D * wxgll(i)
+ weight = jacobian1D * wxgll(i)
elseif(iedge_acoustic ==ILEFT)then
xgamma = - xiz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
zgamma = + xix(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
jacobian1D = sqrt(xgamma**2 + zgamma**2)
nx = - zgamma / jacobian1D
nz = + xgamma / jacobian1D
- weight = jacobian1D * wzgll(j)
+ weight = jacobian1D * wzgll(j)
elseif(iedge_acoustic ==IRIGHT)then
xgamma = - xiz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
zgamma = + xix(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
jacobian1D = sqrt(xgamma**2 + zgamma**2)
nx = + zgamma / jacobian1D
nz = - xgamma / jacobian1D
- weight = jacobian1D * wzgll(j)
+ weight = jacobian1D * wzgll(j)
endif
! compute dot product [u_s + w]*n
@@ -4288,9 +4331,9 @@
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + weight*displ_n
if(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) +&
- weight*((b_displ_x + b_displw_x)*nx + (b_displ_z + b_displw_z)*nz)
- endif !if(SIMULATION_TYPE == 2) then
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
+ + weight*((b_displ_x + b_displw_x)*nx + (b_displ_z + b_displw_z)*nz)
+ endif
enddo
@@ -4303,95 +4346,111 @@
! ************************************ add force source
! ************************************************************************************
- if(any_acoustic) then
+ if(any_acoustic) then
! --- add the source
- if(.not. initialfield) then
+ if(.not. initialfield) then
- do i_source=1,NSOURCES
-! if this processor carries the source and the source element is acoustic
- if (is_proc_source(i_source) == 1 .and. .not. elastic(ispec_selected_source(i_source)) .and. &
- .not. poroelastic(ispec_selected_source(i_source))) then
+ do i_source=1,NSOURCES
+ ! if this processor carries the source and the source element is acoustic
+ if (is_proc_source(i_source) == 1 .and. &
+ .not. elastic(ispec_selected_source(i_source)) .and. &
+ .not. poroelastic(ispec_selected_source(i_source))) then
+
! collocated force
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
! to add minus the source to Chi_dot_dot to get plus the source in pressure
- if(source_type(i_source) == 1) then
+ if(source_type(i_source) == 1) then
- if(SIMULATION_TYPE == 1) then ! forward wavefield
- do j = 1,NGLLZ
- do i = 1,NGLLX
- iglob = ibool(i,j,ispec_selected_source(i_source))
- hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - source_time_function(i_source,it)*hlagrange
- enddo
- enddo
- else ! backward wavefield
- do j = 1,NGLLZ
- do i = 1,NGLLX
- iglob = ibool(i,j,ispec_selected_source(i_source))
- hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
- - source_time_function(i_source,NSTEP-it+1)*hlagrange
- enddo
- enddo
- endif
+ if(SIMULATION_TYPE == 1) then
+ ! forward wavefield
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - source_time_function(i_source,it)*hlagrange
+ enddo
+ enddo
+ else
+ ! backward wavefield
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
+ - source_time_function(i_source,NSTEP-it+1)*hlagrange
+ enddo
+ enddo
+ endif
-! moment tensor
- else if(source_type(i_source) == 2) then
- call exit_MPI('cannot have moment tensor source in acoustic element')
+ ! moment tensor
+ else if(source_type(i_source) == 2) then
+ call exit_MPI('cannot have moment tensor source in acoustic element')
- endif
- endif ! if this processor carries the source and the source element is acoustic
- enddo ! do i_source=1,NSOURCES
+ endif
+ endif ! if this processor carries the source and the source element is acoustic
+ enddo ! do i_source=1,NSOURCES
- if(SIMULATION_TYPE == 2) then ! adjoint wavefield
- irec_local = 0
- do irec = 1,nrec
-! add the source (only if this proc carries the source)
- if (myrank == which_proc_receiver(irec)) then
+ if(SIMULATION_TYPE == 2) then ! adjoint wavefield
+ irec_local = 0
+ do irec = 1,nrec
+ ! add the source (only if this proc carries the source)
+ if (myrank == which_proc_receiver(irec)) then
- irec_local = irec_local + 1
- if (.not. elastic(ispec_selected_rec(irec)) .and. &
- .not. poroelastic(ispec_selected_rec(irec))) then
-! add source array
- do j=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,ispec_selected_rec(irec))
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - &
- adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
- enddo
- enddo
- endif ! if element acoustic
+ irec_local = irec_local + 1
+ if (.not. elastic(ispec_selected_rec(irec)) .and. &
+ .not. poroelastic(ispec_selected_rec(irec))) then
+ ! add source array
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec_selected_rec(irec))
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
+ enddo
+ enddo
+ endif ! if element acoustic
- endif ! if this processor carries the adjoint source
- enddo ! irec = 1,nrec
- endif ! SIMULATION_TYPE == 2 adjoint wavefield
+ endif ! if this processor carries the adjoint source
+ enddo ! irec = 1,nrec
+ endif ! SIMULATION_TYPE == 2 adjoint wavefield
- endif ! if not using an initial field
- endif !if(any_acoustic)
+ endif ! if not using an initial field
+
+ endif !if(any_acoustic)
! assembling potential_dot_dot for acoustic elements
#ifdef USE_MPI
if ( nproc > 1 .and. any_acoustic .and. ninterface_acoustic > 0) then
call assemble_MPI_vector_ac(potential_dot_dot_acoustic,npoin, &
- ninterface, ninterface_acoustic,inum_interfaces_acoustic, &
- max_interface_size, max_ibool_interfaces_size_ac,&
- ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
- tab_requests_send_recv_acoustic,buffer_send_faces_vector_ac, &
- buffer_recv_faces_vector_ac, my_neighbours)
+ ninterface, ninterface_acoustic,inum_interfaces_acoustic, &
+ max_interface_size, max_ibool_interfaces_size_ac,&
+ ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+ tab_requests_send_recv_acoustic,buffer_send_faces_vector_ac, &
+ buffer_recv_faces_vector_ac, my_neighbours)
+
+ if ( SIMULATION_TYPE == 2) then
+ call assemble_MPI_vector_ac(b_potential_dot_dot_acoustic,npoin, &
+ ninterface, ninterface_acoustic,inum_interfaces_acoustic, &
+ max_interface_size, max_ibool_interfaces_size_ac,&
+ ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+ tab_requests_send_recv_acoustic,buffer_send_faces_vector_ac, &
+ buffer_recv_faces_vector_ac, my_neighbours)
+
+ endif
+
endif
- if ( nproc > 1 .and. any_acoustic .and. ninterface_acoustic > 0 .and. SIMULATION_TYPE == 2) then
- call assemble_MPI_vector_ac(b_potential_dot_dot_acoustic,npoin, &
- ninterface, ninterface_acoustic,inum_interfaces_acoustic, &
- max_interface_size, max_ibool_interfaces_size_ac,&
- ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
- tab_requests_send_recv_acoustic,buffer_send_faces_vector_ac, &
- buffer_recv_faces_vector_ac, my_neighbours)
- endif
+! if ( nproc > 1 .and. any_acoustic .and. ninterface_acoustic > 0 .and. SIMULATION_TYPE == 2) then
+! call assemble_MPI_vector_ac(b_potential_dot_dot_acoustic,npoin, &
+! ninterface, ninterface_acoustic,inum_interfaces_acoustic, &
+! max_interface_size, max_ibool_interfaces_size_ac,&
+! ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+! tab_requests_send_recv_acoustic,buffer_send_faces_vector_ac, &
+! buffer_recv_faces_vector_ac, my_neighbours)
+! endif
#endif
! ************************************************************************************
More information about the CIG-COMMITS
mailing list