[cig-commits] r17949 - seismo/2D/SPECFEM2D/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Feb 22 19:14:04 PST 2011
Author: danielpeter
Date: 2011-02-22 19:14:03 -0800 (Tue, 22 Feb 2011)
New Revision: 17949
Modified:
seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
seismo/2D/SPECFEM2D/trunk/check_stability.F90
seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90
seismo/2D/SPECFEM2D/trunk/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_viscoelastic.f90
seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90
seismo/2D/SPECFEM2D/trunk/constants.h.in
seismo/2D/SPECFEM2D/trunk/invert_mass_matrix.f90
seismo/2D/SPECFEM2D/trunk/locate_receivers.F90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
seismo/2D/SPECFEM2D/trunk/prepare_assemble_MPI.F90
seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
updates parallelizations and domain counters in SPECFEM2D/
Modified: seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -54,7 +54,9 @@
!-----------------------------------------------
! Assembling the mass matrix.
!-----------------------------------------------
- subroutine assemble_MPI_scalar(array_val1, array_val2, array_val3, array_val4,npoin, &
+ subroutine assemble_MPI_scalar(array_val1,npoin_val1, &
+ array_val2,npoin_val2, &
+ array_val3,array_val4,npoin_val3, &
ninterface, max_interface_size, max_ibool_interfaces_size_ac, &
max_ibool_interfaces_size_el, &
max_ibool_interfaces_size_po, &
@@ -68,7 +70,6 @@
include 'constants.h'
include 'mpif.h'
- 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, &
@@ -79,7 +80,15 @@
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
+ ! acoustic
+ integer :: npoin_val1
+ real(kind=CUSTOM_REAL), dimension(npoin_val1), intent(inout) :: array_val1
+ ! elastic
+ integer :: npoin_val2
+ real(kind=CUSTOM_REAL), dimension(npoin_val2), intent(inout) :: array_val2
+ ! poroelastic
+ integer :: npoin_val3
+ real(kind=CUSTOM_REAL), dimension(npoin_val3), intent(inout) :: array_val3,array_val4
integer :: ipoin, num_interface
integer :: ier
@@ -88,7 +97,7 @@
2*max_ibool_interfaces_size_po, ninterface) :: &
buffer_send_faces_scalar, &
buffer_recv_faces_scalar
- integer :: msg_status(MPI_STATUS_SIZE)
+ integer, dimension(MPI_STATUS_SIZE) :: msg_status
integer, dimension(ninterface) :: msg_requests
buffer_send_faces_scalar(:,:) = 0.d0
@@ -120,6 +129,7 @@
array_val4(ibool_interfaces_poroelastic(i,num_interface))
end do
+ ! non-blocking synchronous send request
call MPI_ISSEND( buffer_send_faces_scalar(1,num_interface), &
nibool_interfaces_acoustic(num_interface)+nibool_interfaces_elastic(num_interface)+&
nibool_interfaces_poroelastic(num_interface)+nibool_interfaces_poroelastic(num_interface), &
@@ -130,6 +140,8 @@
end do
do num_interface = 1, ninterface
+
+ ! starts a blocking receive
call MPI_recv ( buffer_recv_faces_scalar(1,num_interface), &
nibool_interfaces_acoustic(num_interface)+nibool_interfaces_elastic(num_interface)+&
nibool_interfaces_poroelastic(num_interface)+nibool_interfaces_poroelastic(num_interface), &
@@ -167,6 +179,7 @@
end do
+ ! synchronizes MPI processes
call MPI_BARRIER(mpi_comm_world,ier)
end subroutine assemble_MPI_scalar
@@ -216,45 +229,50 @@
integer, dimension(ninterface), intent(in) :: my_neighbours
! local parameters
- integer :: ipoin, num_interface, inum_interface
- integer :: ier
- integer :: i
+ integer :: ipoin, num_interface,iinterface,ier,iglob
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
+ tab_requests_send_recv_acoustic(:) = 0
- do inum_interface = 1, ninterface_acoustic
+ ! loops over acoustic interfaces only
+ do iinterface = 1, ninterface_acoustic
- num_interface = inum_interfaces_acoustic(inum_interface)
+ ! gets interface index in the range of all interfaces [1,ninterface]
+ num_interface = inum_interfaces_acoustic(iinterface)
- ipoin = 0
- do i = 1, nibool_interfaces_acoustic(num_interface)
- ipoin = ipoin + 1
- buffer_send_faces_vector_ac(ipoin,inum_interface) = &
- array_val1(ibool_interfaces_acoustic(i,num_interface))
- end do
+ ! loops over all interface points
+ do ipoin = 1, nibool_interfaces_acoustic(num_interface)
+ iglob = ibool_interfaces_acoustic(ipoin,num_interface)
+ ! copies array values to buffer
+ buffer_send_faces_vector_ac(ipoin,iinterface) = array_val1(iglob)
+ end do
+
end do
- do inum_interface = 1, ninterface_acoustic
+ do iinterface = 1, ninterface_acoustic
- num_interface = inum_interfaces_acoustic(inum_interface)
+ ! gets global interface index
+ num_interface = inum_interfaces_acoustic(iinterface)
- call MPI_ISSEND( buffer_send_faces_vector_ac(1,inum_interface), &
+ ! non-blocking synchronous send
+ call MPI_ISSEND( buffer_send_faces_vector_ac(1,iinterface), &
nibool_interfaces_acoustic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
- tab_requests_send_recv_acoustic(inum_interface), ier)
+ tab_requests_send_recv_acoustic(iinterface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_ISSEND unsuccessful in assemble_MPI_vector_start')
end if
- call MPI_Irecv ( buffer_recv_faces_vector_ac(1,inum_interface), &
+ ! starts a non-blocking receive
+ call MPI_Irecv ( buffer_recv_faces_vector_ac(1,iinterface), &
nibool_interfaces_acoustic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
- tab_requests_send_recv_acoustic(ninterface_acoustic+inum_interface), ier)
+ tab_requests_send_recv_acoustic(ninterface_acoustic+iinterface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_Irecv unsuccessful in assemble_MPI_vector')
@@ -262,26 +280,37 @@
end do
- do inum_interface = 1, ninterface_acoustic*2
-
- call MPI_Wait (tab_requests_send_recv_acoustic(inum_interface), status_acoustic, ier)
-
+
+ ! waits for MPI requests to complete (recv)
+ ! each wait returns once the specified MPI request completed
+ do iinterface = 1, ninterface_acoustic
+ call MPI_Wait (tab_requests_send_recv_acoustic(ninterface_acoustic+iinterface), &
+ status_acoustic, ier)
enddo
- do inum_interface = 1, ninterface_acoustic
+ ! assembles the array values
+ do iinterface = 1, ninterface_acoustic
- num_interface = inum_interfaces_acoustic(inum_interface)
+ ! gets global interface index
+ num_interface = inum_interfaces_acoustic(iinterface)
- 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)
- end do
+ ! loops over all interface points
+ do ipoin = 1, nibool_interfaces_acoustic(num_interface)
+ iglob = ibool_interfaces_acoustic(ipoin,num_interface)
+ ! adds buffer contribution
+ array_val1(iglob) = array_val1(iglob) + buffer_recv_faces_vector_ac(ipoin,iinterface)
+ end do
end do
+
+ ! waits for MPI requests to complete (send)
+ ! just to make sure that all sending is done
+ do iinterface = 1, ninterface_acoustic
+ call MPI_Wait (tab_requests_send_recv_acoustic(iinterface), status_acoustic, ier)
+ enddo
+
+
end subroutine assemble_MPI_vector_ac
@@ -297,15 +326,14 @@
! communication scheme.
!-----------------------------------------------
subroutine assemble_MPI_vector_el(array_val2,npoin, &
- ninterface, ninterface_elastic, &
- inum_interfaces_elastic, &
- max_interface_size, max_ibool_interfaces_size_el,&
- ibool_interfaces_elastic, nibool_interfaces_elastic, &
- tab_requests_send_recv_elastic, &
- buffer_send_faces_vector_el, &
- buffer_recv_faces_vector_el, &
- my_neighbours &
- )
+ ninterface, ninterface_elastic, &
+ inum_interfaces_elastic, &
+ max_interface_size, max_ibool_interfaces_size_el,&
+ ibool_interfaces_elastic, nibool_interfaces_elastic, &
+ tab_requests_send_recv_elastic, &
+ buffer_send_faces_vector_el, &
+ buffer_recv_faces_vector_el, &
+ my_neighbours)
implicit none
@@ -329,43 +357,40 @@
real(kind=CUSTOM_REAL), dimension(3,npoin), intent(inout) :: array_val2
integer, dimension(ninterface), intent(in) :: my_neighbours
- integer :: ipoin, num_interface, inum_interface
- integer :: ier
+ integer :: ipoin, num_interface, iinterface, ier, i
integer, dimension(MPI_STATUS_SIZE) :: status_elastic
- integer :: i
+ do iinterface = 1, ninterface_elastic
- do inum_interface = 1, ninterface_elastic
+ num_interface = inum_interfaces_elastic(iinterface)
- num_interface = inum_interfaces_elastic(inum_interface)
-
ipoin = 0
do i = 1, nibool_interfaces_elastic(num_interface)
- buffer_send_faces_vector_el(ipoin+1:ipoin+3,inum_interface) = &
+ buffer_send_faces_vector_el(ipoin+1:ipoin+3,iinterface) = &
array_val2(:,ibool_interfaces_elastic(i,num_interface))
ipoin = ipoin + 3
end do
end do
- do inum_interface = 1, ninterface_elastic
+ do iinterface = 1, ninterface_elastic
- num_interface = inum_interfaces_elastic(inum_interface)
+ num_interface = inum_interfaces_elastic(iinterface)
- call MPI_ISSEND( buffer_send_faces_vector_el(1,inum_interface), &
+ call MPI_ISSEND( buffer_send_faces_vector_el(1,iinterface), &
3*nibool_interfaces_elastic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
- tab_requests_send_recv_elastic(inum_interface), ier)
+ tab_requests_send_recv_elastic(iinterface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_ISSEND unsuccessful in assemble_MPI_vector_el')
end if
- call MPI_Irecv ( buffer_recv_faces_vector_el(1,inum_interface), &
+ call MPI_Irecv ( buffer_recv_faces_vector_el(1,iinterface), &
3*nibool_interfaces_elastic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
- tab_requests_send_recv_elastic(ninterface_elastic+inum_interface), ier)
+ tab_requests_send_recv_elastic(ninterface_elastic+iinterface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_Irecv unsuccessful in assemble_MPI_vector_el')
@@ -373,21 +398,21 @@
end do
- do inum_interface = 1, ninterface_elastic*2
+ do iinterface = 1, ninterface_elastic*2
- call MPI_Wait (tab_requests_send_recv_elastic(inum_interface), status_elastic, ier)
+ call MPI_Wait (tab_requests_send_recv_elastic(iinterface), status_elastic, ier)
enddo
- do inum_interface = 1, ninterface_elastic
+ do iinterface = 1, ninterface_elastic
- num_interface = inum_interfaces_elastic(inum_interface)
+ num_interface = inum_interfaces_elastic(iinterface)
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)
+ + buffer_recv_faces_vector_el(ipoin+1:ipoin+3,iinterface)
ipoin = ipoin + 3
end do
@@ -408,15 +433,14 @@
! communication scheme.
!-----------------------------------------------
subroutine assemble_MPI_vector_po(array_val3,array_val4,npoin, &
- ninterface, ninterface_poroelastic, &
- inum_interfaces_poroelastic, &
- max_interface_size, max_ibool_interfaces_size_po,&
- ibool_interfaces_poroelastic, nibool_interfaces_poroelastic, &
- tab_requests_send_recv_poro, &
- buffer_send_faces_vector_pos,buffer_send_faces_vector_pow, &
- buffer_recv_faces_vector_pos,buffer_recv_faces_vector_pow, &
- my_neighbours &
- )
+ ninterface, ninterface_poroelastic, &
+ inum_interfaces_poroelastic, &
+ max_interface_size, max_ibool_interfaces_size_po,&
+ ibool_interfaces_poroelastic, nibool_interfaces_poroelastic, &
+ tab_requests_send_recv_poro, &
+ buffer_send_faces_vector_pos,buffer_send_faces_vector_pow, &
+ buffer_recv_faces_vector_pos,buffer_recv_faces_vector_pow, &
+ my_neighbours)
implicit none
@@ -440,68 +464,65 @@
real(kind=CUSTOM_REAL), dimension(NDIM,npoin), intent(inout) :: array_val3,array_val4
integer, dimension(ninterface), intent(in) :: my_neighbours
- integer :: ipoin, num_interface, inum_interface
- integer :: ier
+ integer :: ipoin, num_interface, iinterface, ier, i
integer, dimension(MPI_STATUS_SIZE) :: status_poroelastic
- integer :: i
+ do iinterface = 1, ninterface_poroelastic
- do inum_interface = 1, ninterface_poroelastic
+ num_interface = inum_interfaces_poroelastic(iinterface)
- num_interface = inum_interfaces_poroelastic(inum_interface)
-
ipoin = 0
do i = 1, nibool_interfaces_poroelastic(num_interface)
- buffer_send_faces_vector_pos(ipoin+1:ipoin+2,inum_interface) = &
+ buffer_send_faces_vector_pos(ipoin+1:ipoin+2,iinterface) = &
array_val3(:,ibool_interfaces_poroelastic(i,num_interface))
ipoin = ipoin + 2
end do
ipoin = 0
do i = 1, nibool_interfaces_poroelastic(num_interface)
- buffer_send_faces_vector_pow(ipoin+1:ipoin+2,inum_interface) = &
+ buffer_send_faces_vector_pow(ipoin+1:ipoin+2,iinterface) = &
array_val4(:,ibool_interfaces_poroelastic(i,num_interface))
ipoin = ipoin + 2
end do
end do
- do inum_interface = 1, ninterface_poroelastic
+ do iinterface = 1, ninterface_poroelastic
- num_interface = inum_interfaces_poroelastic(inum_interface)
+ num_interface = inum_interfaces_poroelastic(iinterface)
- call MPI_ISSEND( buffer_send_faces_vector_pos(1,inum_interface), &
+ call MPI_ISSEND( buffer_send_faces_vector_pos(1,iinterface), &
NDIM*nibool_interfaces_poroelastic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
- tab_requests_send_recv_poro(inum_interface), ier)
+ tab_requests_send_recv_poro(iinterface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_ISSEND unsuccessful in assemble_MPI_vector_pos')
end if
- call MPI_Irecv ( buffer_recv_faces_vector_pos(1,inum_interface), &
+ call MPI_Irecv ( buffer_recv_faces_vector_pos(1,iinterface), &
NDIM*nibool_interfaces_poroelastic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
- tab_requests_send_recv_poro(ninterface_poroelastic+inum_interface), ier)
+ tab_requests_send_recv_poro(ninterface_poroelastic+iinterface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_Irecv unsuccessful in assemble_MPI_vector_pos')
end if
- call MPI_ISSEND( buffer_send_faces_vector_pow(1,inum_interface), &
+ call MPI_ISSEND( buffer_send_faces_vector_pow(1,iinterface), &
NDIM*nibool_interfaces_poroelastic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
- tab_requests_send_recv_poro(ninterface_poroelastic*2+inum_interface), ier)
+ tab_requests_send_recv_poro(ninterface_poroelastic*2+iinterface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_ISSEND unsuccessful in assemble_MPI_vector_pow')
end if
- call MPI_Irecv ( buffer_recv_faces_vector_pow(1,inum_interface), &
+ call MPI_Irecv ( buffer_recv_faces_vector_pow(1,iinterface), &
NDIM*nibool_interfaces_poroelastic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
- tab_requests_send_recv_poro(ninterface_poroelastic*3+inum_interface), ier)
+ tab_requests_send_recv_poro(ninterface_poroelastic*3+iinterface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_Irecv unsuccessful in assemble_MPI_vector_pow')
@@ -509,21 +530,21 @@
end do
- do inum_interface = 1, ninterface_poroelastic*4
+ do iinterface = 1, ninterface_poroelastic*4
- call MPI_Wait (tab_requests_send_recv_poro(inum_interface), status_poroelastic, ier)
+ call MPI_Wait (tab_requests_send_recv_poro(iinterface), status_poroelastic, ier)
enddo
- do inum_interface = 1, ninterface_poroelastic
+ do iinterface = 1, ninterface_poroelastic
- num_interface = inum_interfaces_poroelastic(inum_interface)
+ num_interface = inum_interfaces_poroelastic(iinterface)
ipoin = 0
do i = 1, nibool_interfaces_poroelastic(num_interface)
array_val3(:,ibool_interfaces_poroelastic(i,num_interface)) = &
array_val3(:,ibool_interfaces_poroelastic(i,num_interface)) + &
- buffer_recv_faces_vector_pos(ipoin+1:ipoin+2,inum_interface)
+ buffer_recv_faces_vector_pos(ipoin+1:ipoin+2,iinterface)
ipoin = ipoin + 2
end do
@@ -531,7 +552,7 @@
do i = 1, nibool_interfaces_poroelastic(num_interface)
array_val4(:,ibool_interfaces_poroelastic(i,num_interface)) = &
array_val4(:,ibool_interfaces_poroelastic(i,num_interface)) + &
- buffer_recv_faces_vector_pow(ipoin+1:ipoin+2,inum_interface)
+ buffer_recv_faces_vector_pow(ipoin+1:ipoin+2,iinterface)
ipoin = ipoin + 2
end do
Modified: seismo/2D/SPECFEM2D/trunk/check_stability.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/check_stability.F90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/check_stability.F90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -43,7 +43,8 @@
!========================================================================
- subroutine check_stability(myrank,time,it,NSTEP,npoin, &
+ subroutine check_stability(myrank,time,it,NSTEP, &
+ npoin_acoustic,npoin_elastic,npoin_poroelastic, &
any_elastic_glob,any_elastic,displ_elastic, &
any_poroelastic_glob,any_poroelastic, &
displs_poroelastic,displw_poroelastic, &
@@ -59,18 +60,20 @@
#endif
integer :: myrank,it,NSTEP
- integer :: npoin
double precision :: time
logical :: any_elastic_glob,any_elastic
- real(kind=CUSTOM_REAL), dimension(3,npoin) :: displ_elastic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin_elastic) :: displ_elastic
logical :: any_poroelastic_glob,any_poroelastic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,displw_poroelastic
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin_poroelastic) :: displs_poroelastic,displw_poroelastic
logical :: any_acoustic_glob,any_acoustic
- real(kind=CUSTOM_REAL), dimension(npoin) :: potential_acoustic
+ integer :: npoin_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin_acoustic) :: potential_acoustic
double precision :: time_start
integer :: year_start,month_start
Modified: seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -41,8 +41,10 @@
!
!========================================================================
-subroutine compute_curl_one_element(curl_element,displ_elastic,displs_poroelastic,elastic,poroelastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec)
+ subroutine compute_curl_one_element(curl_element,displ_elastic, &
+ displs_poroelastic,elastic,poroelastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin_elastic,npoin_poroelastic,ispec)
! compute curl in (poro)elastic elements (for rotational study)
@@ -50,7 +52,7 @@
include "constants.h"
- integer nspec,npoin,ispec
+ integer nspec,ispec
integer, dimension(NGLLX,NGLLX,nspec) :: ibool
@@ -60,8 +62,12 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: curl_element
logical, dimension(nspec) :: elastic,poroelastic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displ_elastic,displs_poroelastic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin_elastic) :: displ_elastic
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin_poroelastic) :: displs_poroelastic
+
! array with derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
@@ -92,9 +98,9 @@
! 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)
+ duz_dxi = duz_dxi + displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
- duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
enddo
xixl = xix(i,j,ispec)
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -42,18 +42,19 @@
!
!========================================================================
- subroutine compute_energy(displ_elastic,veloc_elastic,displs_poroelastic,velocs_poroelastic, &
+ 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, &
+ nspec,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ assign_external_model,it,deltat,t0,kmato,poroelastcoef,density, &
porosity,tortuosity, &
vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext, &
anisotropic,anisotropy,wxgll,wzgll,numat, &
pressure_element,vector_field_element,e1,e11, &
potential_dot_acoustic,potential_dot_dot_acoustic, &
- TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS,p_sv)
! compute kinetic and potential energy in the solid (acoustic elements are excluded)
@@ -61,7 +62,7 @@
include "constants.h"
- integer :: nspec,npoin,numat
+ integer :: nspec,numat
! vector field in an element
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
@@ -73,10 +74,11 @@
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) :: &
+ integer :: npoin_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin_acoustic) :: &
potential_dot_acoustic,potential_dot_dot_acoustic
- logical :: TURN_ATTENUATION_ON
+ logical :: TURN_ATTENUATION_ON,p_sv
integer :: it
double precision :: t0,deltat
@@ -86,22 +88,22 @@
logical, dimension(nspec) :: elastic,poroelastic,anisotropic
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
-
integer, dimension(nspec) :: kmato
-
logical :: assign_external_model
-
double precision, dimension(2,numat) :: density
double precision, dimension(numat) :: porosity,tortuosity
double precision, dimension(6,numat) :: anisotropy
- double precision, dimension(4,3,numat) :: elastcoef
+ double precision, dimension(4,3,numat) :: poroelastcoef
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,vsext,rhoext
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
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displw_poroelastic,velocw_poroelastic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin_elastic) :: displ_elastic,veloc_elastic
+
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin_poroelastic) :: displs_poroelastic,velocs_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin_poroelastic) :: displw_poroelastic,velocw_poroelastic
! Gauss-Lobatto-Legendre points and weights
real(kind=CUSTOM_REAL), dimension(NGLLX) :: wxgll
@@ -143,10 +145,15 @@
!---
if(elastic(ispec)) then
+ ! checks wave type
+ if( .not. p_sv ) then
+ call exit_MPI('output energy for SH waves not implemented yet')
+ endif
+
! 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))
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec))
+ lambdalplus2mul_relaxed = poroelastcoef(3,1,kmato(ispec))
rhol = density(1,kmato(ispec))
! double loop over GLL points
@@ -164,19 +171,19 @@
endif
! derivative along x and along z
- dux_dxi = ZERO
- duz_dxi = ZERO
+ dux_dxi = 0._CUSTOM_REAL
+ duz_dxi = 0._CUSTOM_REAL
- dux_dgamma = ZERO
- duz_dgamma = ZERO
+ dux_dgamma = 0._CUSTOM_REAL
+ duz_dgamma = 0._CUSTOM_REAL
! 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)
+ duz_dxi = duz_dxi + displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
- duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
enddo
xixl = xix(i,j,ispec)
@@ -195,7 +202,7 @@
! 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
+ + veloc_elastic(3,ibool(i,j,ispec))**2) *wxgll(i)*wzgll(j)*jacobianl / TWO
! compute potential energy
potential_energy = potential_energy &
@@ -217,15 +224,15 @@
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
+ mul_s = poroelastcoef(2,1,kmato(ispec))
+ kappal_s = poroelastcoef(3,1,kmato(ispec)) - FOUR_THIRDS*mul_s
rhol_s = density(1,kmato(ispec))
!fluid properties
- kappal_f = elastcoef(1,2,kmato(ispec))
+ kappal_f = poroelastcoef(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
+ mul_fr = poroelastcoef(2,3,kmato(ispec))
+ kappal_fr = poroelastcoef(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))
@@ -337,24 +344,25 @@
! 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, &
+ call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic, &
+ displs_poroelastic,displw_poroelastic,elastic,poroelastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec, &
+ npoin_acoustic,npoin_elastic,npoin_poroelastic,assign_external_model, &
+ numat,kmato,density,porosity,tortuosity,poroelastcoef,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)
+ call compute_vector_one_element(vector_field_element,potential_dot_acoustic, &
+ veloc_elastic,velocs_poroelastic, &
+ elastic,poroelastic,xix,xiz,gammax,gammaz, &
+ ibool,hprime_xx,hprime_zz, &
+ nspec,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ ispec,numat,kmato,density,rhoext,assign_external_model)
! get density of current spectral element
- lambdal_relaxed = elastcoef(1,1,kmato(ispec))
- mul_relaxed = elastcoef(2,1,kmato(ispec))
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(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)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_viscoelastic.f90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_viscoelastic.f90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -48,7 +48,7 @@
initialfield,TURN_ATTENUATION_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
accel_elastic,veloc_elastic,displ_elastic,b_accel_elastic,b_displ_elastic, &
- density,elastcoef,xix,xiz,gammax,gammaz, &
+ density,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
source_time_function,sourcearray,adj_sourcearrays,e1,e11, &
e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
@@ -96,7 +96,7 @@
real(kind=CUSTOM_REAL), dimension(3,npoin) :: accel_elastic,veloc_elastic,displ_elastic
double precision, dimension(2,numat) :: density
- double precision, dimension(4,3,numat) :: elastcoef
+ double precision, dimension(4,3,numat) :: poroelastcoef
double precision, dimension(6,numat) :: anisotropy
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,vsext,rhoext
@@ -213,9 +213,9 @@
if(elastic(ispec)) then
! 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))
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec))
+ lambdalplus2mul_relaxed = poroelastcoef(3,1,kmato(ispec))
! first double loop over GLL points to compute and store gradients
do j = 1,NGLLZ
@@ -468,8 +468,8 @@
ispec = numabs(ispecabs)
! get elastic parameters of current spectral element
- lambdal_relaxed = elastcoef(1,1,kmato(ispec))
- mul_relaxed = elastcoef(2,1,kmato(ispec))
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(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)
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -43,11 +43,12 @@
!========================================================================
subroutine compute_pressure_whole_medium(potential_dot_dot_acoustic,displ_elastic,&
- displs_poroelastic,displw_poroelastic,elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext, &
- c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,e1,e11, &
- TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
+ displs_poroelastic,displw_poroelastic,elastic,poroelastic,vector_field_display, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec, &
+ npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic,assign_external_model, &
+ numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext, &
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,e1,e11, &
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute pressure in acoustic elements and in elastic elements
@@ -63,7 +64,7 @@
double precision, dimension(2,numat) :: density
double precision, dimension(numat) :: porosity,tortuosity
- double precision, dimension(4,3,numat) :: elastcoef
+ double precision, dimension(4,3,numat) :: poroelastcoef
double precision, dimension(6,numat) :: anisotropy
double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,vsext,rhoext
double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
@@ -71,9 +72,13 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
logical, dimension(nspec) :: elastic,poroelastic,anisotropic
- real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_dot_acoustic
- real(kind=CUSTOM_REAL), dimension(3,npoin) :: displ_elastic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,displw_poroelastic
+ integer :: npoin_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin_acoustic) :: potential_dot_dot_acoustic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin_elastic) :: displ_elastic
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin_poroelastic) :: displs_poroelastic,displw_poroelastic
+
double precision, dimension(3,npoin) :: vector_field_display
! array with derivatives of Lagrange polynomials
@@ -98,8 +103,9 @@
! compute pressure in this element
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,&
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec, &
+ npoin_acoustic,npoin_elastic,npoin_poroelastic,assign_external_model, &
+ numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext, &
c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,ispec,e1,e11, &
TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
@@ -121,8 +127,9 @@
subroutine compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,&
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec, &
+ npoin_acoustic,npoin_elastic,npoin_poroelastic,assign_external_model, &
+ numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext, &
c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,ispec,e1,e11, &
TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
@@ -132,14 +139,14 @@
include "constants.h"
- integer nspec,npoin,numat,ispec
+ integer nspec,numat,ispec
integer, dimension(nspec) :: kmato
integer, dimension(NGLLX,NGLLX,nspec) :: ibool
double precision, dimension(2,numat) :: density
double precision, dimension(numat) :: porosity,tortuosity
- double precision, dimension(4,3,numat) :: elastcoef
+ double precision, dimension(4,3,numat) :: poroelastcoef
double precision, dimension(6,numat) :: anisotropy
double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,vsext,rhoext
double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
@@ -150,9 +157,12 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: pressure_element
logical, dimension(nspec) :: elastic,poroelastic,anisotropic
- real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_dot_acoustic
- real(kind=CUSTOM_REAL), dimension(3,npoin) :: displ_elastic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,displw_poroelastic
+ integer :: npoin_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin_acoustic) :: potential_dot_dot_acoustic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin_elastic) :: displ_elastic
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin_poroelastic) :: displs_poroelastic,displw_poroelastic
! array with derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
@@ -219,9 +229,9 @@
if(elastic(ispec)) then
! 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))
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec))
+ lambdalplus2mul_relaxed = poroelastcoef(3,1,kmato(ispec))
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -338,22 +348,22 @@
elseif(poroelastic(ispec)) then
- lambdal_relaxed = elastcoef(1,1,kmato(ispec))
- mul_relaxed = elastcoef(2,1,kmato(ispec))
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec))
! get poroelastic 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
+ mul_s = poroelastcoef(2,1,kmato(ispec))
+ kappal_s = poroelastcoef(3,1,kmato(ispec)) - FOUR_THIRDS*mul_s
rhol_s = density(1,kmato(ispec))
!fluid properties
- kappal_f = elastcoef(1,2,kmato(ispec))
+ kappal_f = poroelastcoef(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
+ mul_fr = poroelastcoef(2,3,kmato(ispec))
+ kappal_fr = poroelastcoef(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))
Modified: seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -43,8 +43,10 @@
!========================================================================
subroutine compute_vector_whole_medium(potential_acoustic,veloc_elastic,velocs_poroelastic,&
- elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
+ elastic,poroelastic,vector_field_display, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ numat,kmato,density,rhoext,assign_external_model)
! compute Grad(potential) in acoustic elements
! and combine with existing velocity vector field in elastic elements
@@ -56,21 +58,20 @@
integer nspec,npoin,numat
logical :: assign_external_model
-
integer, dimension(nspec) :: kmato
-
double precision, dimension(NGLLX,NGLLX,nspec) :: rhoext
-
double precision, dimension(2,numat) :: density
-
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
-
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
logical, dimension(nspec) :: elastic,poroelastic
- real(kind=CUSTOM_REAL), dimension(npoin) :: potential_acoustic
- real(kind=CUSTOM_REAL), dimension(3,npoin) :: veloc_elastic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: velocs_poroelastic
+ integer :: npoin_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin_acoustic) :: potential_acoustic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin_elastic) :: veloc_elastic
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin_poroelastic) :: velocs_poroelastic
+
double precision, dimension(3,npoin) :: vector_field_display
! array with derivatives of Lagrange polynomials
@@ -87,9 +88,12 @@
do ispec = 1,nspec
! compute vector field in this element
- call compute_vector_one_element(vector_field_element,potential_acoustic,veloc_elastic,velocs_poroelastic,&
- elastic,poroelastic,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,&
- density,rhoext,assign_external_model)
+ call compute_vector_one_element(vector_field_element,potential_acoustic, &
+ veloc_elastic,velocs_poroelastic, &
+ elastic,poroelastic,xix,xiz,gammax,gammaz, &
+ ibool,hprime_xx,hprime_zz, &
+ nspec,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ ispec,numat,kmato,density,rhoext,assign_external_model)
! store the result
do j = 1,NGLLZ
@@ -107,9 +111,12 @@
!=====================================================================
!
- subroutine compute_vector_one_element(vector_field_element,potential_acoustic,veloc_elastic,velocs_poroelastic,&
- elastic,poroelastic,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,&
- density,rhoext,assign_external_model)
+ subroutine compute_vector_one_element(vector_field_element,potential_acoustic, &
+ veloc_elastic,velocs_poroelastic,&
+ elastic,poroelastic,xix,xiz,gammax,gammaz, &
+ ibool,hprime_xx,hprime_zz, &
+ nspec,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ ispec,numat,kmato,density,rhoext,assign_external_model)
! compute Grad(potential) if acoustic element or copy existing vector if elastic element
@@ -117,7 +124,7 @@
include "constants.h"
- integer nspec,npoin,ispec,numat
+ integer nspec,ispec,numat
logical :: assign_external_model
@@ -135,9 +142,12 @@
real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLX) :: vector_field_element
logical, dimension(nspec) :: elastic,poroelastic
- real(kind=CUSTOM_REAL), dimension(npoin) :: potential_acoustic
- real(kind=CUSTOM_REAL), dimension(3,npoin) :: veloc_elastic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: velocs_poroelastic
+ integer :: npoin_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin_acoustic) :: potential_acoustic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin_elastic) :: veloc_elastic
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin_poroelastic) :: velocs_poroelastic
! array with derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
@@ -173,7 +183,7 @@
do i = 1,NGLLX
iglob = ibool(i,j,ispec)
vector_field_element(1,i,j) = velocs_poroelastic(1,iglob)
- vector_field_element(2,i,j) = ZERO
+ vector_field_element(2,i,j) = 0._CUSTOM_REAL
vector_field_element(3,i,j) = velocs_poroelastic(2,iglob)
enddo
enddo
@@ -189,7 +199,7 @@
do i = 1,NGLLX
! derivative along x
- tempx1l = ZERO
+ tempx1l = 0._CUSTOM_REAL
do k = 1,NGLLX
hp1 = hprime_xx(i,k)
iglob = ibool(k,j,ispec)
@@ -197,7 +207,7 @@
enddo
! derivative along z
- tempx2l = ZERO
+ tempx2l = 0._CUSTOM_REAL
do k = 1,NGLLZ
hp2 = hprime_zz(j,k)
iglob = ibool(i,k,ispec)
@@ -213,7 +223,7 @@
! derivatives of potential
vector_field_element(1,i,j) = (tempx1l*xixl + tempx2l*gammaxl) / rhol
- vector_field_element(2,i,j) = ZERO
+ vector_field_element(2,i,j) = 0._CUSTOM_REAL
vector_field_element(3,i,j) = (tempx1l*xizl + tempx2l*gammazl) / rhol
enddo
Modified: seismo/2D/SPECFEM2D/trunk/constants.h.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h.in 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/constants.h.in 2011-02-23 03:14:03 UTC (rev 17949)
@@ -9,7 +9,7 @@
!
! solver in single or double precision depending on the machine (4 or 8 bytes)
!
-! ALSO CHANGE FILE precision.h ACCORDINGLY
+! ALSO CHANGE FILE precision_mpi.h ACCORDINGLY
!
integer, parameter :: SIZE_REAL = 4
integer, parameter :: SIZE_DOUBLE = 8
@@ -18,7 +18,7 @@
! set to SIZE_REAL to run in single precision
! set to SIZE_DOUBLE to run in double precision (increases memory size by 2)
!
-! DO NOT forget to change precision_mpi.h accordingly
+! DO CHANGE precision_mpi.h accordingly
!
integer, parameter :: CUSTOM_REAL = @CUSTOM_REAL@
Modified: seismo/2D/SPECFEM2D/trunk/invert_mass_matrix.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/invert_mass_matrix.f90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/invert_mass_matrix.f90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -43,11 +43,11 @@
!
!========================================================================
- subroutine invert_mass_matrix_init(any_elastic,any_acoustic,any_poroelastic,npoin, &
- rmass_inverse_elastic,&
- rmass_inverse_acoustic, &
+ subroutine invert_mass_matrix_init(any_elastic,any_acoustic,any_poroelastic, &
+ rmass_inverse_elastic,npoin_elastic, &
+ rmass_inverse_acoustic,npoin_acoustic, &
rmass_s_inverse_poroelastic, &
- rmass_w_inverse_poroelastic, &
+ rmass_w_inverse_poroelastic,npoin_poroelastic, &
nspec,ibool,kmato,wxgll,wzgll,jacobian, &
elastic,poroelastic, &
assign_external_model,numat, &
@@ -60,12 +60,18 @@
include 'constants.h'
logical any_elastic,any_acoustic,any_poroelastic
- integer npoin
! inverse mass matrices
- real(kind=CUSTOM_REAL), dimension(npoin) :: rmass_inverse_elastic,rmass_inverse_acoustic
- real(kind=CUSTOM_REAL), dimension(npoin) :: rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(npoin_elastic) :: rmass_inverse_elastic
+ integer :: npoin_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin_acoustic) :: rmass_inverse_acoustic
+
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(npoin_poroelastic) :: &
+ rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
+
integer :: nspec
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
integer, dimension(nspec) :: kmato
@@ -88,10 +94,10 @@
double precision :: rhol_s,rhol_f,rhol_bar,phil,tortl
! initializes mass matrix
- if(any_elastic) rmass_inverse_elastic(:) = ZERO
- if(any_poroelastic) rmass_s_inverse_poroelastic(:) = ZERO
- if(any_poroelastic) rmass_w_inverse_poroelastic(:) = ZERO
- if(any_acoustic) rmass_inverse_acoustic(:) = ZERO
+ if(any_elastic) rmass_inverse_elastic(:) = 0._CUSTOM_REAL
+ if(any_poroelastic) rmass_s_inverse_poroelastic(:) = 0._CUSTOM_REAL
+ if(any_poroelastic) rmass_w_inverse_poroelastic(:) = 0._CUSTOM_REAL
+ if(any_acoustic) rmass_inverse_acoustic(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do j = 1,NGLLZ
@@ -152,11 +158,11 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic,npoin, &
- rmass_inverse_elastic,&
- rmass_inverse_acoustic, &
+ subroutine invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic, &
+ rmass_inverse_elastic,npoin_elastic, &
+ rmass_inverse_acoustic,npoin_acoustic, &
rmass_s_inverse_poroelastic, &
- rmass_w_inverse_poroelastic)
+ rmass_w_inverse_poroelastic,npoin_poroelastic)
! inverts the global mass matrix
@@ -165,23 +171,36 @@
logical any_elastic,any_acoustic,any_poroelastic
- integer npoin
-
! inverse mass matrices
- real(kind=CUSTOM_REAL), dimension(npoin) :: rmass_inverse_elastic,rmass_inverse_acoustic
- real(kind=CUSTOM_REAL), dimension(npoin) :: rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
+ integer :: npoin_elastic
+ real(kind=CUSTOM_REAL), dimension(npoin_elastic) :: rmass_inverse_elastic
+ integer :: npoin_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin_acoustic) :: rmass_inverse_acoustic
+
+ integer :: npoin_poroelastic
+ real(kind=CUSTOM_REAL), dimension(npoin_poroelastic) :: &
+ rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
+
! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
- if(any_elastic) where(rmass_inverse_elastic <= 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
- if(any_poroelastic) where(rmass_s_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_s_inverse_poroelastic = 1._CUSTOM_REAL
- if(any_poroelastic) where(rmass_w_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_w_inverse_poroelastic = 1._CUSTOM_REAL
- if(any_acoustic) where(rmass_inverse_acoustic <= 0._CUSTOM_REAL) rmass_inverse_acoustic = 1._CUSTOM_REAL
+ if(any_elastic) &
+ where(rmass_inverse_elastic <= 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
+ if(any_poroelastic) &
+ where(rmass_s_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_s_inverse_poroelastic = 1._CUSTOM_REAL
+ if(any_poroelastic) &
+ where(rmass_w_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_w_inverse_poroelastic = 1._CUSTOM_REAL
+ if(any_acoustic) &
+ where(rmass_inverse_acoustic <= 0._CUSTOM_REAL) rmass_inverse_acoustic = 1._CUSTOM_REAL
! compute the inverse of the mass matrix
- if(any_elastic) rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
- if(any_poroelastic) rmass_s_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_s_inverse_poroelastic(:)
- if(any_poroelastic) rmass_w_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_w_inverse_poroelastic(:)
- if(any_acoustic) rmass_inverse_acoustic(:) = 1._CUSTOM_REAL / rmass_inverse_acoustic(:)
+ if(any_elastic) &
+ rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
+ if(any_poroelastic) &
+ rmass_s_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_s_inverse_poroelastic(:)
+ if(any_poroelastic) &
+ rmass_w_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_w_inverse_poroelastic(:)
+ if(any_acoustic) &
+ rmass_inverse_acoustic(:) = 1._CUSTOM_REAL / rmass_inverse_acoustic(:)
end subroutine invert_mass_matrix
Modified: seismo/2D/SPECFEM2D/trunk/locate_receivers.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_receivers.F90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/locate_receivers.F90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -46,10 +46,13 @@
!---- locate_receivers finds the correct position of the receivers
!----
- subroutine locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank, &
- st_xval,st_zval,ispec_selected_rec, &
- xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo,ipass, &
- x_final_receiver, z_final_receiver)
+ subroutine locate_receivers(ibool,coord,nspec,npoin,xigll,zigll, &
+ nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank, &
+ st_xval,st_zval,ispec_selected_rec, &
+ xi_receiver,gamma_receiver,station_name,network_name, &
+ x_source,z_source, &
+ coorg,knods,ngnod,npgeo,ipass, &
+ x_final_receiver, z_final_receiver)
implicit none
@@ -107,6 +110,9 @@
ierror = 0
+#ifdef USE_MPI
+ call MPI_BARRIER(MPI_COMM_WORLD,ierror)
+#endif
! **************
@@ -128,177 +134,179 @@
! loop on all the stations
do irec=1,nrec
-! set distance to huge initial value
- distmin=HUGEVAL
+ ! set distance to huge initial value
+ distmin=HUGEVAL
read(1,*) station_name(irec),network_name(irec),st_xval(irec),st_zval(irec),stele,stbur
-! check that station is not buried, burial is not implemented in current code
+ ! check that station is not buried, burial is not implemented in current code
if(abs(stbur) > TINYVAL) call exit_MPI('stations with non-zero burial not implemented yet')
-! compute distance between source and receiver
- distance_receiver(irec) = sqrt((st_zval(irec)-z_source)**2 + (st_xval(irec)-x_source)**2)
+ ! compute distance between source and receiver
+ distance_receiver(irec) = sqrt((st_zval(irec)-z_source)**2 + (st_xval(irec)-x_source)**2)
- do ispec=1,nspec
+ do ispec=1,nspec
-! loop only on points inside the element
-! exclude edges to ensure this point is not shared with other elements
- do j=2,NGLLZ-1
- do i=2,NGLLX-1
+ ! loop only on points inside the element
+ ! exclude edges to ensure this point is not shared with other elements
+ do j=2,NGLLZ-1
+ do i=2,NGLLX-1
- iglob = ibool(i,j,ispec)
- dist = sqrt((st_xval(irec)-dble(coord(1,iglob)))**2 + (st_zval(irec)-dble(coord(2,iglob)))**2)
+ iglob = ibool(i,j,ispec)
+ dist = sqrt((st_xval(irec)-dble(coord(1,iglob)))**2 + (st_zval(irec)-dble(coord(2,iglob)))**2)
-! keep this point if it is closer to the receiver
- if(dist < distmin) then
- distmin = dist
- ispec_selected_rec(irec) = ispec
- ix_initial_guess = i
- iz_initial_guess = j
- endif
+ ! keep this point if it is closer to the receiver
+ if(dist < distmin) then
+ distmin = dist
+ ispec_selected_rec(irec) = ispec
+ ix_initial_guess = i
+ iz_initial_guess = j
+ endif
- enddo
enddo
-
-! end of loop on all the spectral elements
enddo
+ ! end of loop on all the spectral elements
+ enddo
+
! ****************************************
! find the best (xi,gamma) for each receiver
! ****************************************
-! use initial guess in xi and gamma
- xi = xigll(ix_initial_guess)
- gamma = zigll(iz_initial_guess)
+ ! use initial guess in xi and gamma
+ xi = xigll(ix_initial_guess)
+ gamma = zigll(iz_initial_guess)
-! iterate to solve the non linear system
- do iter_loop = 1,NUM_ITER
+ ! iterate to solve the non linear system
+ do iter_loop = 1,NUM_ITER
-! recompute jacobian for the new point
- call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_rec(irec),ngnod,nspec,npgeo, &
- .true.)
+ ! recompute jacobian for the new point
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian, &
+ coorg,knods,ispec_selected_rec(irec),ngnod,nspec,npgeo, &
+ .true.)
-! compute distance to target location
- dx = - (x - st_xval(irec))
- dz = - (z - st_zval(irec))
+ ! compute distance to target location
+ dx = - (x - st_xval(irec))
+ dz = - (z - st_zval(irec))
-! compute increments
- dxi = xix*dx + xiz*dz
- dgamma = gammax*dx + gammaz*dz
+ ! compute increments
+ dxi = xix*dx + xiz*dz
+ dgamma = gammax*dx + gammaz*dz
-! update values
- xi = xi + dxi
- gamma = gamma + dgamma
+ ! update values
+ xi = xi + dxi
+ gamma = gamma + dgamma
-! impose that we stay in that element
-! (useful if user gives a receiver outside the mesh for instance)
-! we can go slightly outside the [1,1] segment since with finite elements
-! the polynomial solution is defined everywhere
-! this can be useful for convergence of itertive scheme with distorted elements
- if (xi > 1.10d0) xi = 1.10d0
- if (xi < -1.10d0) xi = -1.10d0
- if (gamma > 1.10d0) gamma = 1.10d0
- if (gamma < -1.10d0) gamma = -1.10d0
+ ! impose that we stay in that element
+ ! (useful if user gives a receiver outside the mesh for instance)
+ ! we can go slightly outside the [1,1] segment since with finite elements
+ ! the polynomial solution is defined everywhere
+ ! this can be useful for convergence of itertive scheme with distorted elements
+ if (xi > 1.10d0) xi = 1.10d0
+ if (xi < -1.10d0) xi = -1.10d0
+ if (gamma > 1.10d0) gamma = 1.10d0
+ if (gamma < -1.10d0) gamma = -1.10d0
-! end of non linear iterations
- enddo
+ ! end of non linear iterations
+ enddo
-! compute final coordinates of point found
- call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_rec(irec),ngnod,nspec,npgeo, &
- .true.)
+ ! compute final coordinates of point found
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian, &
+ coorg,knods,ispec_selected_rec(irec),ngnod,nspec,npgeo, &
+ .true.)
-! store xi,gamma of point found
- xi_receiver(irec) = xi
- gamma_receiver(irec) = gamma
+ ! store xi,gamma of point found
+ xi_receiver(irec) = xi
+ gamma_receiver(irec) = gamma
-! compute final distance between asked and found
- final_distance(irec) = sqrt((st_xval(irec)-x)**2 + (st_zval(irec)-z)**2)
+ ! compute final distance between asked and found
+ final_distance(irec) = sqrt((st_xval(irec)-x)**2 + (st_zval(irec)-z)**2)
- x_final_receiver(irec) = x
- z_final_receiver(irec) = z
+ x_final_receiver(irec) = x
+ z_final_receiver(irec) = z
-enddo
+ enddo
-! close receiver file
-close(1)
+ ! close receiver file
+ close(1)
! elect one process for each receiver.
#ifdef USE_MPI
-call MPI_GATHER(final_distance(1),nrec,MPI_DOUBLE_PRECISION,&
- gather_final_distance(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
-call MPI_GATHER(xi_receiver(1),nrec,MPI_DOUBLE_PRECISION,&
- gather_xi_receiver(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
-call MPI_GATHER(gamma_receiver(1),nrec,MPI_DOUBLE_PRECISION,&
- gather_gamma_receiver(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
-call MPI_GATHER(ispec_selected_rec(1),nrec,MPI_INTEGER,&
- gather_ispec_selected_rec(1,1),nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
+ call MPI_GATHER(final_distance(1),nrec,MPI_DOUBLE_PRECISION,&
+ gather_final_distance(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
+ call MPI_GATHER(xi_receiver(1),nrec,MPI_DOUBLE_PRECISION,&
+ gather_xi_receiver(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
+ call MPI_GATHER(gamma_receiver(1),nrec,MPI_DOUBLE_PRECISION,&
+ gather_gamma_receiver(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
+ call MPI_GATHER(ispec_selected_rec(1),nrec,MPI_INTEGER,&
+ gather_ispec_selected_rec(1,1),nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
-if ( myrank == 0 ) then
- do irec = 1, nrec
+ if ( myrank == 0 ) then
+ do irec = 1, nrec
which_proc_receiver(irec:irec) = minloc(gather_final_distance(irec,:)) - 1
+ enddo
+ endif
- enddo
-endif
+ call MPI_BCAST(which_proc_receiver(1),nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
-call MPI_BCAST(which_proc_receiver(1),nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
-
#else
-gather_final_distance(:,1) = final_distance(:)
+ gather_final_distance(:,1) = final_distance(:)
-gather_xi_receiver(:,1) = xi_receiver(:)
-gather_gamma_receiver(:,1) = gamma_receiver(:)
-gather_ispec_selected_rec(:,1) = ispec_selected_rec(:)
+ gather_xi_receiver(:,1) = xi_receiver(:)
+ gather_gamma_receiver(:,1) = gamma_receiver(:)
+ gather_ispec_selected_rec(:,1) = ispec_selected_rec(:)
-which_proc_receiver(:) = 0
+ which_proc_receiver(:) = 0
#endif
-nrecloc = 0
-do irec = 1, nrec
- if ( which_proc_receiver(irec) == myrank ) then
+ nrecloc = 0
+ do irec = 1, nrec
+ if ( which_proc_receiver(irec) == myrank ) then
nrecloc = nrecloc + 1
recloc(nrecloc) = irec
- endif
-enddo
+ endif
+ enddo
-if (myrank == 0 .and. ipass == 1) then
+ if (myrank == 0 .and. ipass == 1) then
- do irec = 1, nrec
- write(IOUT,*)
- write(IOUT,*) 'Station # ',irec,' ',station_name(irec),network_name(irec)
+ do irec = 1, nrec
+ write(IOUT,*)
+ write(IOUT,*) 'Station # ',irec,' ',station_name(irec),network_name(irec)
- if(gather_final_distance(irec,which_proc_receiver(irec)+1) == HUGEVAL) call exit_MPI('error locating receiver')
+ if(gather_final_distance(irec,which_proc_receiver(irec)+1) == HUGEVAL) &
+ call exit_MPI('error locating receiver')
- write(IOUT,*) ' original x: ',sngl(st_xval(irec))
- write(IOUT,*) ' original z: ',sngl(st_zval(irec))
- write(IOUT,*) ' distance from source: ',sngl(distance_receiver(irec))
- write(IOUT,*) 'closest estimate found: ',sngl(gather_final_distance(irec,which_proc_receiver(irec)+1)),' m away'
- write(IOUT,*) ' in element ',gather_ispec_selected_rec(irec,which_proc_receiver(irec)+1)
- write(IOUT,*) ' at process ', which_proc_receiver(irec)
- write(IOUT,*) ' at xi,gamma coordinates = ',gather_xi_receiver(irec,which_proc_receiver(irec)+1),&
- gather_gamma_receiver(irec,which_proc_receiver(irec)+1)
+ write(IOUT,*) ' original x: ',sngl(st_xval(irec))
+ write(IOUT,*) ' original z: ',sngl(st_zval(irec))
+ write(IOUT,*) ' distance from source: ',sngl(distance_receiver(irec))
+ write(IOUT,*) 'closest estimate found: ',sngl(gather_final_distance(irec,which_proc_receiver(irec)+1)), &
+ ' m away'
+ write(IOUT,*) ' in element ',gather_ispec_selected_rec(irec,which_proc_receiver(irec)+1)
+ write(IOUT,*) ' at process ', which_proc_receiver(irec)
+ write(IOUT,*) ' at xi,gamma coordinates = ',gather_xi_receiver(irec,which_proc_receiver(irec)+1),&
+ gather_gamma_receiver(irec,which_proc_receiver(irec)+1)
+ write(IOUT,*)
+ enddo
+
write(IOUT,*)
+ write(IOUT,*) 'end of receiver detection'
+ write(IOUT,*)
- enddo
-
- write(IOUT,*)
- write(IOUT,*) 'end of receiver detection'
- write(IOUT,*)
-
- ! write out actual station locations (compare with STATIONS_target from meshfem2D)
- ! NOTE: this will be written out even if generate_STATIONS = .false.
- open(unit=15,file='DATA/STATIONS',status='unknown')
- do irec = 1,nrec
- write(15,"('S',i4.4,' AA ',f20.7,1x,f20.7,' 0.0 0.0')") &
+ ! write out actual station locations (compare with STATIONS_target from meshfem2D)
+ ! NOTE: this will be written out even if generate_STATIONS = .false.
+ open(unit=15,file='DATA/STATIONS',status='unknown')
+ do irec = 1,nrec
+ write(15,"('S',i4.4,' AA ',f20.7,1x,f20.7,' 0.0 0.0')") &
irec,x_final_receiver(irec),z_final_receiver(irec)
- enddo
- close(15)
+ enddo
+ close(15)
-endif
+ endif
-! deallocate arrays
+ ! deallocate arrays
deallocate(final_distance)
#ifdef USE_MPI
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -758,11 +758,20 @@
!*****************************
! partitioning
!*****************************
-
+
+ ! allocates & initializes partioning of elements
allocate(part(0:nelmnts-1))
+ part(:) = -1
+
+ if( nproc > 1 ) then
+ allocate(xadj_g(0:nelmnts))
+ allocate(adjncy_g(0:MAX_NEIGHBORS*nelmnts-1))
+ xadj_g(:) = 0
+ adjncy_g(:) = -1
+ endif
! construction of the graph
-
+
! if ngnod == 9, we work on a subarray of elements that represents the elements with four nodes (four corners) only
! because the adjacency of the mesh elements can be entirely determined from the knowledge of the four corners only
if ( ngnod == 9 ) then
@@ -777,22 +786,25 @@
!! DK DK (nxread+1)*(nzread+1) is OK for a regular internal mesh only, not for non structured external meshes
!! DK DK call mesh2dual_ncommonnodes(nelmnts, (nxread+1)*(nzread+1), elmnts_bis, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
!! DK DK the subset of element corners is not renumbered therefore we must still use the nnodes computed for 9 nodes here
- call mesh2dual_ncommonnodes(elmnts_bis,1)
+ ! determines maximum neighbors based on 1 common node
+ call mesh2dual_ncommonnodes(elmnts_bis,1,xadj_g,adjncy_g)
endif
else
if ( nproc > 1 ) then
- call mesh2dual_ncommonnodes(elmnts,1)
+ ! determines maximum neighbors based on 1 common node
+ call mesh2dual_ncommonnodes(elmnts,1,xadj_g,adjncy_g)
endif
endif
if ( nproc == 1 ) then
- part(:) = 0
+ part(:) = 0 ! single process has rank 0
else
- nb_edges = xadj(nelmnts)
+ ! number of common edges
+ nb_edges = xadj_g(nelmnts)
! giving weight to edges and vertices. Currently not used.
call read_weights()
@@ -855,33 +867,27 @@
call Construct_glob2loc_elmnts(nproc)
if ( ngnod == 9 ) then
- if ( nproc > 1 ) then
- deallocate(nnodes_elmnts)
- deallocate(nodes_elmnts)
- endif
- allocate(nnodes_elmnts(0:nnodes-1))
- allocate(nodes_elmnts(0:nsize*nnodes-1))
- nnodes_elmnts(:) = 0
- nodes_elmnts(:) = 0
- do i = 0, ngnod*nelmnts-1
+ if( allocated(nnodes_elmnts) ) deallocate(nnodes_elmnts)
+ if( allocated(nodes_elmnts) ) deallocate(nodes_elmnts)
+ allocate(nnodes_elmnts(0:nnodes-1))
+ allocate(nodes_elmnts(0:nsize*nnodes-1))
+ nnodes_elmnts(:) = 0
+ nodes_elmnts(:) = 0
+ do i = 0, ngnod*nelmnts-1
+ nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/ngnod
+ nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
+ enddo
+ else
+ if ( nproc < 2 ) then
+ if( .not. allocated(nnodes_elmnts) ) allocate(nnodes_elmnts(0:nnodes-1))
+ if( .not. allocated(nodes_elmnts) ) allocate(nodes_elmnts(0:nsize*nnodes-1))
+ nnodes_elmnts(:) = 0
+ nodes_elmnts(:) = 0
+ do i = 0, ngnod*nelmnts-1
nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/ngnod
nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
-
- enddo
- else
- if ( nproc < 2 ) then
- allocate(nnodes_elmnts(0:nnodes-1))
- allocate(nodes_elmnts(0:nsize*nnodes-1))
- nnodes_elmnts(:) = 0
- nodes_elmnts(:) = 0
- do i = 0, ngnod*nelmnts-1
- nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/ngnod
- nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
-
- enddo
-
- endif
-
+ enddo
+ endif
endif
! local number of each node for each partition
Modified: seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -53,41 +53,41 @@
integer :: nelmnts
integer, dimension(:), pointer :: elmnts
- integer, dimension(:), pointer :: elmnts_bis
- integer, dimension(:), pointer :: xadj
- integer, dimension(:), pointer :: vwgt
- integer, dimension(:), pointer :: glob2loc_elmnts
- integer, dimension(:), pointer :: part
+ integer, dimension(:), allocatable :: elmnts_bis
+ integer, dimension(:), allocatable :: vwgt
+ integer, dimension(:), allocatable :: glob2loc_elmnts
+ integer, dimension(:), allocatable :: part
integer :: nb_edges
- integer, dimension(:), pointer :: adjwgt
+ integer, dimension(:), allocatable :: adjwgt
- integer, dimension(:), pointer :: adjncy
+ integer, dimension(:), allocatable :: xadj_g
+ integer, dimension(:), allocatable :: adjncy_g
integer :: nnodes
- double precision, dimension(:,:), pointer :: nodes_coords
- integer, dimension(:), pointer :: nnodes_elmnts
- integer, dimension(:), pointer :: nodes_elmnts
- integer, dimension(:), pointer :: glob2loc_nodes_nparts
- integer, dimension(:), pointer :: glob2loc_nodes_parts
- integer, dimension(:), pointer :: glob2loc_nodes
+ double precision, dimension(:,:), allocatable :: nodes_coords
+ integer, dimension(:), allocatable :: nnodes_elmnts
+ integer, dimension(:), allocatable :: nodes_elmnts
+ integer, dimension(:), allocatable :: glob2loc_nodes_nparts
+ integer, dimension(:), allocatable :: glob2loc_nodes_parts
+ integer, dimension(:), allocatable :: glob2loc_nodes
! interface data
integer :: ninterfaces
- integer, dimension(:), pointer :: tab_size_interfaces, tab_interfaces
+ integer, dimension(:), allocatable :: tab_size_interfaces, tab_interfaces
integer :: nelem_acoustic_surface
integer, dimension(:,:), pointer :: acoustic_surface
integer :: nelem_acoustic_surface_loc
integer :: nelemabs
- integer, dimension(:,:), pointer :: abs_surface
- logical, dimension(:,:), pointer :: abs_surface_char
- integer, dimension(:), pointer :: abs_surface_merge
+ integer, dimension(:,:), allocatable :: abs_surface
+ logical, dimension(:,:), allocatable :: abs_surface_char
+ integer, dimension(:), allocatable :: abs_surface_merge
integer :: nelemabs_loc
integer :: nelemabs_merge
- integer, dimension(:), pointer :: ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ integer, dimension(:), allocatable :: ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right
! for acoustic/elastic coupled elements
@@ -326,88 +326,102 @@
!-----------------------------------------------
! Creating dual graph (adjacency is defined by 'ncommonnodes' between two elements).
!-----------------------------------------------
- subroutine mesh2dual_ncommonnodes(elmnts_l,ncommonnodes)
+ subroutine mesh2dual_ncommonnodes(elmnts_l,ncommonnodes,xadj,adjncy)
implicit none
include "constants.h"
integer, dimension(0:NCORNERS*nelmnts-1), intent(in) :: elmnts_l
integer, intent(in) :: ncommonnodes
+ integer, dimension(0:nelmnts),intent(out) :: xadj
+ integer, dimension(0:MAX_NEIGHBORS*nelmnts-1),intent(out) :: adjncy
- integer :: i, j, k, l, m, nb_edges
+ ! local parameters
+ integer :: i, j, k, l, m, num_edges
logical :: is_neighbour
integer :: num_node, n
integer :: elem_base, elem_target
integer :: connectivity
- allocate(xadj(0:nelmnts))
+ ! allocates memory for arrays
+ if( .not. allocated(nnodes_elmnts) ) allocate(nnodes_elmnts(0:nnodes-1))
+ if( .not. allocated(nodes_elmnts) ) allocate(nodes_elmnts(0:nsize*nnodes-1))
+
+ ! initializes
xadj(:) = 0
- allocate(adjncy(0:MAX_NEIGHBORS*nelmnts-1))
adjncy(:) = 0
- allocate(nnodes_elmnts(0:nnodes-1))
nnodes_elmnts(:) = 0
- allocate(nodes_elmnts(0:nsize*nnodes-1))
nodes_elmnts(:) = 0
+ num_edges = 0
- nb_edges = 0
-
! list of elements per node
do i = 0, NCORNERS*nelmnts-1
- nodes_elmnts(elmnts_l(i)*nsize+nnodes_elmnts(elmnts_l(i))) = i/NCORNERS
+ nodes_elmnts(elmnts_l(i)*nsize + nnodes_elmnts(elmnts_l(i))) = i/NCORNERS
nnodes_elmnts(elmnts_l(i)) = nnodes_elmnts(elmnts_l(i)) + 1
enddo
! checking which elements are neighbours ('ncommonnodes' criteria)
do j = 0, nnodes-1
- do k = 0, nnodes_elmnts(j)-1
- do l = k+1, nnodes_elmnts(j)-1
+ do k = 0, nnodes_elmnts(j)-1
+ do l = k+1, nnodes_elmnts(j)-1
- connectivity = 0
- elem_base = nodes_elmnts(k+j*nsize)
- elem_target = nodes_elmnts(l+j*nsize)
- do n = 1, NCORNERS
- num_node = elmnts_l(NCORNERS*elem_base+n-1)
- do m = 0, nnodes_elmnts(num_node)-1
- if ( nodes_elmnts(m+num_node*nsize) == elem_target ) then
- connectivity = connectivity + 1
- endif
- enddo
- enddo
+ connectivity = 0
+ elem_base = nodes_elmnts(k+j*nsize)
+ elem_target = nodes_elmnts(l+j*nsize)
+ do n = 1, NCORNERS
+ num_node = elmnts_l(NCORNERS*elem_base+n-1)
+ do m = 0, nnodes_elmnts(num_node)-1
+ if ( nodes_elmnts(m+num_node*nsize) == elem_target ) then
+ connectivity = connectivity + 1
+ endif
+ enddo
+ enddo
- if ( connectivity >= ncommonnodes) then
+ ! sets adjacency (adjncy) and number of neighbors (xadj)
+ ! according to ncommonnodes criteria
+ if ( connectivity >= ncommonnodes) then
- is_neighbour = .false.
+ is_neighbour = .false.
- do m = 0, xadj(nodes_elmnts(k+j*nsize))
- if ( .not.is_neighbour ) then
- if ( adjncy(nodes_elmnts(k+j*nsize)*MAX_NEIGHBORS+m) == nodes_elmnts(l+j*nsize) ) then
- is_neighbour = .true.
+ do m = 0, xadj(nodes_elmnts(k+j*nsize))
+ if ( .not.is_neighbour ) then
+ if ( adjncy(nodes_elmnts(k+j*nsize)*MAX_NEIGHBORS+m) == nodes_elmnts(l+j*nsize) ) then
+ is_neighbour = .true.
+ endif
+ endif
+ enddo
+ if ( .not.is_neighbour ) then
+ adjncy(nodes_elmnts(k+j*nsize)*MAX_NEIGHBORS &
+ + xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
- endif
- endif
- enddo
- if ( .not.is_neighbour ) then
- adjncy(nodes_elmnts(k+j*nsize)*MAX_NEIGHBORS+xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
- xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
- adjncy(nodes_elmnts(l+j*nsize)*MAX_NEIGHBORS+xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
- xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
- endif
- endif
- enddo
- enddo
+ xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
+ if (xadj(nodes_elmnts(k+j*nsize)) > MAX_NEIGHBORS) &
+ stop 'ERROR : too much neighbours per element, modify the mesh.'
+
+ adjncy(nodes_elmnts(l+j*nsize)*MAX_NEIGHBORS &
+ + xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
+
+ xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
+ if (xadj(nodes_elmnts(l+j*nsize))>MAX_NEIGHBORS) &
+ stop 'ERROR : too much neighbours per element, modify the mesh.'
+
+ endif
+ endif
+ enddo
+ enddo
enddo
! making adjacency arrays compact (to be used for partitioning)
do i = 0, nelmnts-1
- k = xadj(i)
- xadj(i) = nb_edges
- do j = 0, k-1
- adjncy(nb_edges) = adjncy(i*MAX_NEIGHBORS+j)
- nb_edges = nb_edges + 1
- enddo
+ k = xadj(i)
+ xadj(i) = num_edges
+ do j = 0, k-1
+ adjncy(num_edges) = adjncy(i*MAX_NEIGHBORS+j)
+ num_edges = num_edges + 1
+ enddo
enddo
- xadj(nelmnts) = nb_edges
+ xadj(nelmnts) = num_edges
end subroutine mesh2dual_ncommonnodes
@@ -442,16 +456,16 @@
allocate(glob2loc_elmnts(0:nelmnts-1))
+ ! initializes number of local elements per partition
do num_part = 0, nparts-1
- num_loc(num_part) = 0
-
+ num_loc(num_part) = 0
enddo
+ ! local numbering
do num_glob = 0, nelmnts-1
- num_part = part(num_glob)
- glob2loc_elmnts(num_glob) = num_loc(num_part)
- num_loc(num_part) = num_loc(num_part) + 1
-
+ num_part = part(num_glob)
+ glob2loc_elmnts(num_glob) = num_loc(num_part)
+ num_loc(num_part) = num_loc(num_part) + 1
enddo
end subroutine Construct_glob2loc_elmnts
@@ -485,14 +499,12 @@
glob2loc_nodes_nparts(num_node) = size_glob2loc_nodes
do el = 0, nnodes_elmnts(num_node)-1
parts_node(part(nodes_elmnts(el+nsize*num_node))) = 1
-
enddo
do num_part = 0, nparts-1
if ( parts_node(num_part) == 1 ) then
size_glob2loc_nodes = size_glob2loc_nodes + 1
parts_node(num_part) = 0
-
endif
enddo
@@ -513,7 +525,6 @@
do num_node = 0, nnodes-1
do el = 0, nnodes_elmnts(num_node)-1
parts_node(part(nodes_elmnts(el+nsize*num_node))) = 1
-
enddo
do num_part = 0, nparts-1
@@ -572,42 +583,53 @@
do num_part_bis = num_part+1, nparts-1
do el = 0, nelmnts-1
if ( part(el) == num_part ) then
+ ! sets material flag
if ( phi_material(num_material(el+1)) < TINYVAL) then
- is_acoustic_el = .false.
- is_elastic_el = .true.
+ ! elastic element
+ is_acoustic_el = .false.
+ is_elastic_el = .true.
elseif ( phi_material(num_material(el+1)) >= 1.d0) then
- is_acoustic_el = .true.
- is_elastic_el = .false.
+ ! acoustic element
+ is_acoustic_el = .true.
+ is_elastic_el = .false.
else
- is_acoustic_el = .false.
- is_elastic_el = .false.
+ ! poroelastic element
+ is_acoustic_el = .false.
+ is_elastic_el = .false.
endif
- do el_adj = xadj(el), xadj(el+1)-1
- if ( phi_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
- is_acoustic_el_adj = .false.
- is_elastic_el_adj = .true.
- elseif ( phi_material(num_material(adjncy(el_adj)+1)) >= 1.d0) then
- is_acoustic_el_adj = .true.
- is_elastic_el_adj = .false.
- else
- is_acoustic_el_adj = .false.
- is_elastic_el_adj = .false.
- endif
- if ( (part(adjncy(el_adj)) == num_part_bis) .and. (is_acoustic_el .eqv. is_acoustic_el_adj) &
- .and. (is_elastic_el .eqv. is_elastic_el_adj) ) then
+
+ ! looks at all neighbor elements
+ do el_adj = xadj_g(el), xadj_g(el+1)-1
+ ! sets neighbor material flag
+ if ( phi_material(num_material(adjncy_g(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .false.
+ is_elastic_el_adj = .true.
+ elseif ( phi_material(num_material(adjncy_g(el_adj)+1)) >= 1.d0) then
+ is_acoustic_el_adj = .true.
+ is_elastic_el_adj = .false.
+ else
+ is_acoustic_el_adj = .false.
+ is_elastic_el_adj = .false.
+ endif
+ ! adds element if neighbor element lies in next parition
+ ! and belongs to same material
+ if ( (part(adjncy_g(el_adj)) == num_part_bis) .and. &
+ (is_acoustic_el .eqv. is_acoustic_el_adj) .and. &
+ (is_elastic_el .eqv. is_elastic_el_adj) ) then
num_edge = num_edge + 1
-
- endif
+ endif
enddo
endif
enddo
+ ! stores number of elements at interface
tab_size_interfaces(num_interface+1) = tab_size_interfaces(num_interface) + num_edge
num_edge = 0
num_interface = num_interface + 1
enddo
enddo
-
+
+ ! stores element indices for elements from above search at each interface
num_interface = 0
num_edge = 0
@@ -615,59 +637,61 @@
tab_interfaces(:) = 0
do num_part = 0, nparts-1
- do num_part_bis = num_part+1, nparts-1
- do el = 0, nelmnts-1
- if ( part(el) == num_part ) then
- if ( phi_material(num_material(el+1)) < TINYVAL) then
- is_acoustic_el = .false.
- is_elastic_el = .true.
- elseif ( phi_material(num_material(el+1)) >= 1.d0) then
- is_acoustic_el = .true.
- is_elastic_el = .false.
+ do num_part_bis = num_part+1, nparts-1
+ do el = 0, nelmnts-1
+ if ( part(el) == num_part ) then
+ if ( phi_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .false.
+ is_elastic_el = .true.
+ elseif ( phi_material(num_material(el+1)) >= 1.d0) then
+ is_acoustic_el = .true.
+ is_elastic_el = .false.
+ else
+ is_acoustic_el = .false.
+ is_elastic_el = .false.
+ endif
+ do el_adj = xadj_g(el), xadj_g(el+1)-1
+ if ( phi_material(num_material(adjncy_g(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .false.
+ is_elastic_el_adj = .true.
+ elseif ( phi_material(num_material(adjncy_g(el_adj)+1)) >= 1.d0) then
+ is_acoustic_el_adj = .true.
+ is_elastic_el_adj = .false.
+ else
+ is_acoustic_el_adj = .false.
+ is_elastic_el_adj = .false.
+ endif
+ if ( (part(adjncy_g(el_adj)) == num_part_bis) .and. &
+ (is_acoustic_el .eqv. is_acoustic_el_adj) .and. &
+ (is_elastic_el .eqv. is_elastic_el_adj) ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+0) = el
+ tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+1) = adjncy_g(el_adj)
+ ncommon_nodes = 0
+ do num_node = 0, 4-1
+ do num_node_bis = 0, 4-1
+ if ( elmnts_l(el*NCORNERS+num_node) == &
+ elmnts_l(adjncy_g(el_adj)*NCORNERS+num_node_bis) ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+3+ncommon_nodes) &
+ = elmnts_l(el*NCORNERS+num_node)
+ ncommon_nodes = ncommon_nodes + 1
+ endif
+ enddo
+ enddo
+ if ( ncommon_nodes > 0 ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+2) = ncommon_nodes
else
- is_acoustic_el = .false.
- is_elastic_el = .false.
+ print *, "Error while building interfaces!", ncommon_nodes
+ stop 'fatal error'
endif
- do el_adj = xadj(el), xadj(el+1)-1
- if ( phi_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
- is_acoustic_el_adj = .false.
- is_elastic_el_adj = .true.
- elseif ( phi_material(num_material(adjncy(el_adj)+1)) >= 1.d0) then
- is_acoustic_el_adj = .true.
- is_elastic_el_adj = .false.
- else
- is_acoustic_el_adj = .false.
- is_elastic_el_adj = .false.
- endif
- if ( (part(adjncy(el_adj)) == num_part_bis) .and. (is_acoustic_el .eqv. is_acoustic_el_adj) &
- .and. (is_elastic_el .eqv. is_elastic_el_adj) ) then
- tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+0) = el
- tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+1) = adjncy(el_adj)
- ncommon_nodes = 0
- do num_node = 0, 4-1
- do num_node_bis = 0, 4-1
- if ( elmnts_l(el*NCORNERS+num_node) == elmnts_l(adjncy(el_adj)*NCORNERS+num_node_bis) ) then
- tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+3+ncommon_nodes) &
- = elmnts_l(el*NCORNERS+num_node)
- ncommon_nodes = ncommon_nodes + 1
- endif
- enddo
- enddo
- if ( ncommon_nodes > 0 ) then
- tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+2) = ncommon_nodes
- else
- print *, "Error while building interfaces!", ncommon_nodes
- stop 'fatal error'
- endif
- num_edge = num_edge + 1
- endif
- enddo
- endif
+ num_edge = num_edge + 1
+ endif
+ enddo
+ endif
- enddo
- num_edge = 0
- num_interface = num_interface + 1
- enddo
+ enddo
+ num_edge = 0
+ num_interface = num_interface + 1
+ enddo
enddo
end subroutine Construct_interfaces
@@ -693,9 +717,7 @@
do j = glob2loc_nodes_nparts(i), glob2loc_nodes_nparts(i+1)-1
if ( glob2loc_nodes_parts(j) == iproc ) then
npgeo = npgeo + 1
-
endif
-
enddo
enddo
else
@@ -1257,6 +1279,7 @@
edgecut = vwgt(0)
edgecut = 0
+ ! we use default strategy for partitioning, thus omit specifing explicit strategy .
call scotchfstratinit (SCOTCHSTRAT(1), IERR)
IF (IERR .NE. 0) THEN
PRINT *, 'ERROR : MAIN : Cannot initialize strat'
@@ -1269,11 +1292,17 @@
STOP
ENDIF
+ ! fills graph structure : see user manual (scotch_user5.1.pdf, page 72/73)
+ ! arguments: #(1) graph_structure #(2) baseval(either 0/1) #(3) number_of_vertices
+ ! #(4) adjacency_index_array #(5) adjacency_end_index_array (optional)
+ ! #(6) vertex_load_array (optional) #(7) vertex_label_array
+ ! #(7) number_of_arcs #(8) adjacency_array
+ ! #(9) arc_load_array (optional) #(10) ierror
CALL SCOTCHFGRAPHBUILD (SCOTCHGRAPH (1), 0, nelmnts, &
- xadj (0), xadj (0), &
- xadj (0), xadj (0), &
- nb_edges, &
- adjncy (0), adjwgt (0), IERR)
+ xadj_g(0), xadj_g(0), &
+ xadj_g(0), xadj_g(0), &
+ nb_edges, &
+ adjncy_g(0), adjwgt (0), IERR)
IF (IERR .NE. 0) THEN
PRINT *, 'ERROR : MAIN : Cannot build graph'
STOP
@@ -1322,15 +1351,20 @@
double precision, dimension(nb_materials), intent(in) :: phi_material
integer, dimension(1:nelmnts), intent(in) :: num_material
-
+ ! local parameters
+ integer, dimension(:), allocatable :: xadj_l
+ integer, dimension(:), allocatable :: adjncy_l
logical, dimension(nb_materials) :: is_acoustic, is_elastic
-
integer :: i, num_edge
integer :: el, el_adj
logical :: is_repartitioned
+ allocate(xadj_l(0:nelmnts))
+ allocate(adjncy_l(0:MAX_NEIGHBORS*nelmnts-1))
+
is_acoustic(:) = .false.
is_elastic(:) = .false.
+
do i = 1, nb_materials
if (phi_material(i) >= 1.d0) then
is_acoustic(i) = .true.
@@ -1340,16 +1374,16 @@
endif
enddo
- call mesh2dual_ncommonnodes(elmnts_l, 2)
+ ! determines maximum neighbors based on 2 common nodes (common edge)
+ call mesh2dual_ncommonnodes(elmnts_l, 2, xadj_l, adjncy_l)
nedges_coupled = 0
do el = 0, nelmnts-1
if ( is_acoustic(num_material(el+1)) ) then
- do el_adj = xadj(el), xadj(el+1) - 1
- if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
+ do el_adj = xadj_l(el), xadj_l(el+1) - 1
+ if ( is_elastic(num_material(adjncy_l(el_adj)+1)) ) then
nedges_coupled = nedges_coupled + 1
endif
-
enddo
endif
enddo
@@ -1359,11 +1393,11 @@
nedges_coupled = 0
do el = 0, nelmnts-1
if ( is_acoustic(num_material(el+1)) ) then
- do el_adj = xadj(el), xadj(el+1) - 1
- if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
+ do el_adj = xadj_l(el), xadj_l(el+1) - 1
+ if ( is_elastic(num_material(adjncy_l(el_adj)+1)) ) then
nedges_coupled = nedges_coupled + 1
edges_coupled(1,nedges_coupled) = el
- edges_coupled(2,nedges_coupled) = adjncy(el_adj)
+ edges_coupled(2,nedges_coupled) = adjncy_l(el_adj)
endif
enddo
@@ -1388,6 +1422,8 @@
endif
enddo
+ deallocate(xadj_l,adjncy_l)
+
end subroutine acoustic_elastic_repartitioning
@@ -1406,15 +1442,20 @@
double precision, dimension(nb_materials), intent(in) :: phi_material
integer, dimension(1:nelmnts), intent(in) :: num_material
-
+ ! local parameters
+ integer, dimension(:), allocatable :: xadj_l
+ integer, dimension(:), allocatable :: adjncy_l
logical, dimension(nb_materials) :: is_acoustic,is_poroelastic
-
integer :: i, num_edge
integer :: el, el_adj
logical :: is_repartitioned
+ allocate(xadj_l(0:nelmnts))
+ allocate(adjncy_l(0:MAX_NEIGHBORS*nelmnts-1))
+
is_acoustic(:) = .false.
is_poroelastic(:) = .false.
+
do i = 1, nb_materials
if (phi_material(i) >=1.d0) then
is_acoustic(i) = .true.
@@ -1424,13 +1465,14 @@
endif
enddo
- call mesh2dual_ncommonnodes(elmnts_l, 2)
+ ! determines maximum neighbors based on 2 common nodes (common edge)
+ call mesh2dual_ncommonnodes(elmnts_l, 2, xadj_l, adjncy_l)
nedges_acporo_coupled = 0
do el = 0, nelmnts-1
if ( is_acoustic(num_material(el+1)) ) then
- do el_adj = xadj(el), xadj(el+1) - 1
- if ( is_poroelastic(num_material(adjncy(el_adj)+1)) ) then
+ do el_adj = xadj_l(el), xadj_l(el+1) - 1
+ if ( is_poroelastic(num_material(adjncy_l(el_adj)+1)) ) then
nedges_acporo_coupled = nedges_acporo_coupled + 1
endif
@@ -1445,11 +1487,11 @@
nedges_acporo_coupled = 0
do el = 0, nelmnts-1
if ( is_acoustic(num_material(el+1)) ) then
- do el_adj = xadj(el), xadj(el+1) - 1
- if ( is_poroelastic(num_material(adjncy(el_adj)+1)) ) then
+ do el_adj = xadj_l(el), xadj_l(el+1) - 1
+ if ( is_poroelastic(num_material(adjncy_l(el_adj)+1)) ) then
nedges_acporo_coupled = nedges_acporo_coupled + 1
edges_acporo_coupled(1,nedges_acporo_coupled) = el
- edges_acporo_coupled(2,nedges_acporo_coupled) = adjncy(el_adj)
+ edges_acporo_coupled(2,nedges_acporo_coupled) = adjncy_l(el_adj)
endif
enddo
@@ -1474,6 +1516,8 @@
endif
enddo
+ deallocate(xadj_l,adjncy_l)
+
end subroutine acoustic_poro_repartitioning
@@ -1492,14 +1536,20 @@
double precision, dimension(nb_materials), intent(in) :: phi_material
integer, dimension(1:nelmnts), intent(in) :: num_material
+ ! local parameters
+ integer, dimension(:), allocatable :: xadj_l
+ integer, dimension(:), allocatable :: adjncy_l
logical, dimension(nb_materials) :: is_elastic,is_poroelastic
-
integer :: i, num_edge
integer :: el, el_adj
logical :: is_repartitioned
+ allocate(xadj_l(0:nelmnts))
+ allocate(adjncy_l(0:MAX_NEIGHBORS*nelmnts-1))
+
is_elastic(:) = .false.
is_poroelastic(:) = .false.
+
do i = 1, nb_materials
if (phi_material(i) < TINYVAL) then
is_elastic(i) = .true.
@@ -1509,13 +1559,14 @@
endif
enddo
- call mesh2dual_ncommonnodes(elmnts_l, 2)
+ ! determines maximum neighbors based on 2 common nodes (common edge)
+ call mesh2dual_ncommonnodes(elmnts_l, 2, xadj_l, adjncy_l)
nedges_elporo_coupled = 0
do el = 0, nelmnts-1
if ( is_poroelastic(num_material(el+1)) ) then
- do el_adj = xadj(el), xadj(el+1) - 1
- if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
+ do el_adj = xadj_l(el), xadj_l(el+1) - 1
+ if ( is_elastic(num_material(adjncy_l(el_adj)+1)) ) then
nedges_elporo_coupled = nedges_elporo_coupled + 1
endif
@@ -1530,11 +1581,11 @@
nedges_elporo_coupled = 0
do el = 0, nelmnts-1
if ( is_poroelastic(num_material(el+1)) ) then
- do el_adj = xadj(el), xadj(el+1) - 1
- if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
+ do el_adj = xadj_l(el), xadj_l(el+1) - 1
+ if ( is_elastic(num_material(adjncy_l(el_adj)+1)) ) then
nedges_elporo_coupled = nedges_elporo_coupled + 1
edges_elporo_coupled(1,nedges_elporo_coupled) = el
- edges_elporo_coupled(2,nedges_elporo_coupled) = adjncy(el_adj)
+ edges_elporo_coupled(2,nedges_elporo_coupled) = adjncy_l(el_adj)
endif
enddo
@@ -1559,6 +1610,8 @@
endif
enddo
+ deallocate(xadj_l,adjncy_l)
+
end subroutine poro_elastic_repartitioning
Modified: seismo/2D/SPECFEM2D/trunk/prepare_assemble_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/prepare_assemble_MPI.F90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/prepare_assemble_MPI.F90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -116,6 +116,7 @@
nibool_interfaces_poroelastic(:) = 0
do num_interface = 1, ninterface
+ ! initializes interface point counters
npoin_interface_acoustic = 0
npoin_interface_elastic = 0
npoin_interface_poroelastic = 0
@@ -197,8 +198,9 @@
ninterface_elastic = 0
ninterface_poroelastic = 0
+ ! loops over all MPI interfaces
do num_interface = 1, ninterface
- ! acoustic
+ ! sets acoustic MPI interface (local) indices in range [1,ninterface_acoustic]
if ( nibool_interfaces_acoustic(num_interface) > 0 ) then
ninterface_acoustic = ninterface_acoustic + 1
inum_interfaces_acoustic(ninterface_acoustic) = num_interface
Modified: seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -155,16 +155,16 @@
endif
-
- ! locate receivers in the mesh
- call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank,&
- st_xval,st_zval,ispec_selected_rec, &
- xi_receiver,gamma_receiver,station_name,network_name,x_source(i_source),z_source(i_source),&
- coorg,knods,ngnod,npgeo,ipass, &
- x_final_receiver, z_final_receiver)
-
enddo ! do i_source=1,NSOURCES
+ ! locate receivers in the mesh
+ call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll, &
+ nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank, &
+ st_xval,st_zval,ispec_selected_rec, &
+ xi_receiver,gamma_receiver,station_name,network_name, &
+ x_source(1),z_source(1), &
+ coorg,knods,ngnod,npgeo,ipass, &
+ x_final_receiver,z_final_receiver)
end subroutine setup_sources_receivers
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -360,6 +360,9 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: curl_element
integer :: i,j,k,it,irec,id,n,ispec,npoin,npgeo,iglob
+ integer :: npoin_acoustic
+ integer :: npoin_elastic
+ integer :: npoin_poroelastic
logical :: anyabs
double precision :: dxd,dyd,dzd,dcurld,valux,valuy,valuz,valcurl,hlagrange,rhol,xi,gamma,x,z
@@ -410,7 +413,8 @@
potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
! inverse mass matrices
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic,rmass_inverse_acoustic
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_acoustic
real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
@@ -1107,38 +1111,7 @@
if( anyabs ) then
-
-! 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
-
-! Files to save absorbed waves needed to reconstruct backward wavefield for adjoint method
+ ! Files to save absorbed waves needed to reconstruct backward wavefield for adjoint method
if(ipass == 1) then
if(any_elastic .and. (SAVE_FORWARD .or. SIMULATION_TYPE == 2)) then
allocate(b_absorb_elastic_left(3,NGLLZ,nspec_xmin,NSTEP))
@@ -1183,23 +1156,8 @@
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
-
else
-! if(.not. allocated(ib_left)) then
-! allocate(ib_left(1))
-! allocate(ib_right(1))
-! allocate(ib_bottom(1))
-! allocate(ib_top(1))
-! endif
-
if(.not. allocated(b_absorb_elastic_left)) then
allocate(b_absorb_elastic_left(1,1,1,1))
allocate(b_absorb_elastic_right(1,1,1,1))
@@ -1464,7 +1422,7 @@
!
!---- build the global mass matrix and invert it once and for all
!
- rmass_inverse_elastic(:) = ZERO
+ rmass_inverse_elastic(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -1492,7 +1450,7 @@
! set entries that are equal to zero to something else, e.g. 1, to avoid division by zero
! these degrees of freedom correspond to points that have been replaced with their periodic counterpart
! and thus are not used any more
- where(rmass_inverse_elastic == ZERO) rmass_inverse_elastic = 1._CUSTOM_REAL
+ where(rmass_inverse_elastic == 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
endif ! of if(ADD_PERIODIC_CONDITIONS)
@@ -2196,17 +2154,16 @@
if(ipass == 1) then
if(any_elastic) then
- allocate(displ_elastic(3,npoin))
- allocate(veloc_elastic(3,npoin))
- allocate(accel_elastic(3,npoin))
- allocate(rmass_inverse_elastic(npoin))
+ npoin_elastic = npoin
else
- ! allocate unused arrays with fictitious size
- allocate(displ_elastic(1,1))
- allocate(veloc_elastic(1,1))
- allocate(accel_elastic(1,1))
- allocate(rmass_inverse_elastic(1))
+ ! allocate unused arrays with fictitious size
+ npoin_elastic = 1
endif
+ allocate(displ_elastic(3,npoin_elastic))
+ allocate(veloc_elastic(3,npoin_elastic))
+ allocate(accel_elastic(3,npoin_elastic))
+ allocate(rmass_inverse_elastic(npoin_elastic))
+
! extra array if adjoint and kernels calculation
if(SIMULATION_TYPE == 2 .and. any_elastic) then
allocate(b_displ_elastic(3,npoin))
@@ -2251,25 +2208,20 @@
endif
if(any_poroelastic) then
- allocate(displs_poroelastic(NDIM,npoin))
- allocate(velocs_poroelastic(NDIM,npoin))
- allocate(accels_poroelastic(NDIM,npoin))
- allocate(rmass_s_inverse_poroelastic(npoin))
- allocate(displw_poroelastic(NDIM,npoin))
- allocate(velocw_poroelastic(NDIM,npoin))
- allocate(accelw_poroelastic(NDIM,npoin))
- allocate(rmass_w_inverse_poroelastic(npoin))
+ npoin_poroelastic = npoin
else
- ! allocate unused arrays with fictitious size
- allocate(displs_poroelastic(1,1))
- allocate(velocs_poroelastic(1,1))
- allocate(accels_poroelastic(1,1))
- allocate(rmass_s_inverse_poroelastic(1))
- allocate(displw_poroelastic(1,1))
- allocate(velocw_poroelastic(1,1))
- allocate(accelw_poroelastic(1,1))
- allocate(rmass_w_inverse_poroelastic(1))
+ ! allocate unused arrays with fictitious size
+ npoin_poroelastic = 1
endif
+ allocate(displs_poroelastic(NDIM,npoin_poroelastic))
+ allocate(velocs_poroelastic(NDIM,npoin_poroelastic))
+ allocate(accels_poroelastic(NDIM,npoin_poroelastic))
+ allocate(rmass_s_inverse_poroelastic(npoin_poroelastic))
+ allocate(displw_poroelastic(NDIM,npoin_poroelastic))
+ allocate(velocw_poroelastic(NDIM,npoin_poroelastic))
+ allocate(accelw_poroelastic(NDIM,npoin_poroelastic))
+ allocate(rmass_w_inverse_poroelastic(npoin_poroelastic))
+
! extra array if adjoint and kernels calculation
if(SIMULATION_TYPE == 2 .and. any_poroelastic) then
allocate(b_displs_poroelastic(NDIM,npoin))
@@ -2375,17 +2327,16 @@
! potential, its first and second derivative, and inverse of the mass matrix for acoustic elements
if(any_acoustic) then
- allocate(potential_acoustic(npoin))
- allocate(potential_dot_acoustic(npoin))
- allocate(potential_dot_dot_acoustic(npoin))
- allocate(rmass_inverse_acoustic(npoin))
+ npoin_acoustic = npoin
else
- ! allocate unused arrays with fictitious size
- allocate(potential_acoustic(1))
- allocate(potential_dot_acoustic(1))
- allocate(potential_dot_dot_acoustic(1))
- allocate(rmass_inverse_acoustic(1))
+ ! allocate unused arrays with fictitious size
+ npoin_acoustic = 1
endif
+ allocate(potential_acoustic(npoin_acoustic))
+ allocate(potential_dot_acoustic(npoin_acoustic))
+ allocate(potential_dot_dot_acoustic(npoin_acoustic))
+ allocate(rmass_inverse_acoustic(npoin_acoustic))
+
if(SIMULATION_TYPE == 2 .and. any_acoustic) then
allocate(b_potential_acoustic(npoin))
allocate(b_potential_dot_acoustic(npoin))
@@ -2424,11 +2375,11 @@
!
!---- build the global mass matrix
!
- call invert_mass_matrix_init(any_elastic,any_acoustic,any_poroelastic,npoin, &
- rmass_inverse_elastic,&
- rmass_inverse_acoustic, &
+ call invert_mass_matrix_init(any_elastic,any_acoustic,any_poroelastic, &
+ rmass_inverse_elastic,npoin_elastic, &
+ rmass_inverse_acoustic,npoin_acoustic, &
rmass_s_inverse_poroelastic, &
- rmass_w_inverse_poroelastic, &
+ rmass_w_inverse_poroelastic,npoin_poroelastic, &
nspec,ibool,kmato,wxgll,wzgll,jacobian, &
elastic,poroelastic, &
assign_external_model,numat, &
@@ -2500,10 +2451,9 @@
endif
! assembling the mass matrix
- call assemble_MPI_scalar(rmass_inverse_acoustic, &
- rmass_inverse_elastic, &
- rmass_s_inverse_poroelastic, &
- rmass_w_inverse_poroelastic,npoin, &
+ call assemble_MPI_scalar(rmass_inverse_acoustic,npoin_acoustic, &
+ rmass_inverse_elastic,npoin_elastic, &
+ rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic,npoin_poroelastic, &
ninterface, max_interface_size, max_ibool_interfaces_size_ac, &
max_ibool_interfaces_size_el, &
max_ibool_interfaces_size_po, &
@@ -2711,8 +2661,11 @@
!--- end of section performed in two passes
!---
- call invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic,npoin,rmass_inverse_elastic,&
- rmass_inverse_acoustic,rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic)
+ call invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic,&
+ rmass_inverse_elastic,npoin_elastic, &
+ rmass_inverse_acoustic,npoin_acoustic, &
+ rmass_s_inverse_poroelastic, &
+ rmass_w_inverse_poroelastic,npoin_poroelastic)
! check the mesh, stability and number of points per wavelength
if(DISPLAY_SUBSET_OPTION == 1) then
@@ -2821,24 +2774,24 @@
!
!---- initialize seismograms
!
- sisux = ZERO
+ sisux = ZERO ! double precision zero
sisuz = ZERO
! initialize arrays to zero
- displ_elastic = ZERO
- veloc_elastic = ZERO
- accel_elastic = ZERO
+ displ_elastic = 0._CUSTOM_REAL
+ veloc_elastic = 0._CUSTOM_REAL
+ accel_elastic = 0._CUSTOM_REAL
- displs_poroelastic = ZERO
- velocs_poroelastic = ZERO
- accels_poroelastic = ZERO
- displw_poroelastic = ZERO
- velocw_poroelastic = ZERO
- accelw_poroelastic = ZERO
+ displs_poroelastic = 0._CUSTOM_REAL
+ velocs_poroelastic = 0._CUSTOM_REAL
+ accels_poroelastic = 0._CUSTOM_REAL
+ displw_poroelastic = 0._CUSTOM_REAL
+ velocw_poroelastic = 0._CUSTOM_REAL
+ accelw_poroelastic = 0._CUSTOM_REAL
- potential_acoustic = ZERO
- potential_dot_acoustic = ZERO
- potential_dot_dot_acoustic = ZERO
+ potential_acoustic = 0._CUSTOM_REAL
+ potential_dot_acoustic = 0._CUSTOM_REAL
+ potential_dot_dot_acoustic = 0._CUSTOM_REAL
!
!----- Files where viscous damping are saved during forward wavefield calculation
@@ -2924,17 +2877,17 @@
open(unit = 98, file = 'OUTPUT_FILES/'//outputname,status = 'unknown',iostat=ios)
if (ios /= 0) stop 'Error writing snapshot to disk'
- rho_kl(:,:,:) = ZERO
- mu_kl(:,:,:) = ZERO
- kappa_kl(:,:,:) = ZERO
+ rho_kl(:,:,:) = 0._CUSTOM_REAL
+ mu_kl(:,:,:) = 0._CUSTOM_REAL
+ kappa_kl(:,:,:) = 0._CUSTOM_REAL
- rhop_kl(:,:,:) = ZERO
- beta_kl(:,:,:) = ZERO
- alpha_kl(:,:,:) = ZERO
- rhorho_el_hessian_final2(:,:,:) = ZERO
- rhorho_el_hessian_temp2(:) = ZERO
- rhorho_el_hessian_final1(:,:,:) = ZERO
- rhorho_el_hessian_temp1(:) = ZERO
+ rhop_kl(:,:,:) = 0._CUSTOM_REAL
+ beta_kl(:,:,:) = 0._CUSTOM_REAL
+ alpha_kl(:,:,:) = 0._CUSTOM_REAL
+ rhorho_el_hessian_final2(:,:,:) = 0._CUSTOM_REAL
+ rhorho_el_hessian_temp2(:) = 0._CUSTOM_REAL
+ rhorho_el_hessian_final1(:,:,:) = 0._CUSTOM_REAL
+ rhorho_el_hessian_temp1(:) = 0._CUSTOM_REAL
endif
if(any_poroelastic) then
@@ -2970,30 +2923,30 @@
open(unit = 19, file = 'OUTPUT_FILES/'//outputname,status = 'unknown',iostat=ios)
if (ios /= 0) stop 'Error writing snapshot to disk'
- rhot_kl(:,:,:) = ZERO
- rhof_kl(:,:,:) = ZERO
- eta_kl(:,:,:) = ZERO
- sm_kl(:,:,:) = ZERO
- mufr_kl(:,:,:) = ZERO
- B_kl(:,:,:) = ZERO
- C_kl(:,:,:) = ZERO
- M_kl(:,:,:) = ZERO
+ rhot_kl(:,:,:) = 0._CUSTOM_REAL
+ rhof_kl(:,:,:) = 0._CUSTOM_REAL
+ eta_kl(:,:,:) = 0._CUSTOM_REAL
+ sm_kl(:,:,:) = 0._CUSTOM_REAL
+ mufr_kl(:,:,:) = 0._CUSTOM_REAL
+ B_kl(:,:,:) = 0._CUSTOM_REAL
+ C_kl(:,:,:) = 0._CUSTOM_REAL
+ M_kl(:,:,:) = 0._CUSTOM_REAL
- rhob_kl(:,:,:) = ZERO
- rhofb_kl(:,:,:) = ZERO
- phi_kl(:,:,:) = ZERO
- mufrb_kl(:,:,:) = ZERO
- Bb_kl(:,:,:) = ZERO
- Cb_kl(:,:,:) = ZERO
- Mb_kl(:,:,:) = ZERO
+ rhob_kl(:,:,:) = 0._CUSTOM_REAL
+ rhofb_kl(:,:,:) = 0._CUSTOM_REAL
+ phi_kl(:,:,:) = 0._CUSTOM_REAL
+ mufrb_kl(:,:,:) = 0._CUSTOM_REAL
+ Bb_kl(:,:,:) = 0._CUSTOM_REAL
+ Cb_kl(:,:,:) = 0._CUSTOM_REAL
+ Mb_kl(:,:,:) = 0._CUSTOM_REAL
- rhobb_kl(:,:,:) = ZERO
- rhofbb_kl(:,:,:) = ZERO
- phib_kl(:,:,:) = ZERO
- cs_kl(:,:,:) = ZERO
- cpI_kl(:,:,:) = ZERO
- cpII_kl(:,:,:) = ZERO
- ratio_kl(:,:,:) = ZERO
+ rhobb_kl(:,:,:) = 0._CUSTOM_REAL
+ rhofbb_kl(:,:,:) = 0._CUSTOM_REAL
+ phib_kl(:,:,:) = 0._CUSTOM_REAL
+ cs_kl(:,:,:) = 0._CUSTOM_REAL
+ cpI_kl(:,:,:) = 0._CUSTOM_REAL
+ cpII_kl(:,:,:) = 0._CUSTOM_REAL
+ ratio_kl(:,:,:) = 0._CUSTOM_REAL
endif
if(any_acoustic) then
@@ -3004,13 +2957,13 @@
open(unit = 96, file = 'OUTPUT_FILES/'//outputname,status = 'unknown',iostat=ios)
if (ios /= 0) stop 'Error writing snapshot to disk'
- rho_ac_kl(:,:,:) = ZERO
- kappa_ac_kl(:,:,:) = ZERO
+ rho_ac_kl(:,:,:) = 0._CUSTOM_REAL
+ kappa_ac_kl(:,:,:) = 0._CUSTOM_REAL
- rhop_ac_kl(:,:,:) = ZERO
- alpha_ac_kl(:,:,:) = ZERO
- rhorho_ac_hessian_final2(:,:,:) = ZERO
- rhorho_ac_hessian_final1(:,:,:) = ZERO
+ rhop_ac_kl(:,:,:) = 0._CUSTOM_REAL
+ alpha_ac_kl(:,:,:) = 0._CUSTOM_REAL
+ rhorho_ac_hessian_final2(:,:,:) = 0._CUSTOM_REAL
+ rhorho_ac_hessian_final1(:,:,:) = 0._CUSTOM_REAL
endif
endif ! if(SIMULATION_TYPE == 2)
@@ -3025,12 +2978,14 @@
if(initialfield) then
! Calculation of the initial field for a plane wave
- call prepare_initialfield(myrank,any_acoustic,any_poroelastic,over_critical_angle, &
+ if( any_elastic ) then
+ call prepare_initialfield(myrank,any_acoustic,any_poroelastic,over_critical_angle, &
NSOURCES,source_type,angleforce,x_source,z_source,f0, &
npoin,numat,poroelastcoef,density,coord, &
angleforce_refl,c_inc,c_refl,cploc,csloc,time_offset, &
A_plane, B_plane, C_plane, &
accel_elastic,veloc_elastic,displ_elastic)
+ endif
if( over_critical_angle ) then
@@ -3108,7 +3063,7 @@
if(.not. initialfield) then
allocate(source_time_function(NSOURCES,NSTEP))
- source_time_function(:,:) = 0.0_CUSTOM_REAL
+ source_time_function(:,:) = 0._CUSTOM_REAL
! computes source time function array
call prepare_source_time_function(myrank,NSTEP,NSOURCES,source_time_function, &
@@ -4501,12 +4456,14 @@
source_time_function,sourcearray,adj_sourcearrays, &
e1,e11,e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
+ hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1, &
+ phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
deltat,coord,add_Bielak_conditions, x0_source, z0_source, &
A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset, f0(1),&
v0x_left(1,it),v0z_left(1,it),v0x_right(1,it),v0z_right(1,it),v0x_bot(1,it),v0z_bot(1,it), &
t0x_left(1,it),t0z_left(1,it),t0x_right(1,it),t0z_right(1,it),t0x_bot(1,it),t0z_bot(1,it), &
- count_left,count_right,count_bottom,over_critical_angle,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
+ count_left,count_right,count_bottom,over_critical_angle, &
+ NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
b_absorb_elastic_left,b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top, &
nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k)
@@ -5848,34 +5805,34 @@
! elastic medium
if(any_elastic) then
- write(outputname,'(a,i6.6,a)') 'lastframe_elastic',myrank,'.bin'
- open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
+ write(outputname,'(a,i6.6,a)') 'lastframe_elastic',myrank,'.bin'
+ open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
if(p_sv)then !P-SV waves
- do j=1,npoin
- read(55) (b_displ_elastic(i,j), i=1,NDIM), &
- (b_veloc_elastic(i,j), i=1,NDIM), &
- (b_accel_elastic(i,j), i=1,NDIM)
- enddo
- b_displ_elastic(3,:) = b_displ_elastic(2,:)
- b_displ_elastic(2,:) = ZERO
- b_veloc_elastic(3,:) = b_veloc_elastic(2,:)
- b_veloc_elastic(2,:) = ZERO
- b_accel_elastic(3,:) = b_accel_elastic(2,:)
- b_accel_elastic(2,:) = ZERO
+ do j=1,npoin
+ read(55) (b_displ_elastic(i,j), i=1,NDIM), &
+ (b_veloc_elastic(i,j), i=1,NDIM), &
+ (b_accel_elastic(i,j), i=1,NDIM)
+ enddo
+ b_displ_elastic(3,:) = b_displ_elastic(2,:)
+ b_displ_elastic(2,:) = 0._CUSTOM_REAL
+ b_veloc_elastic(3,:) = b_veloc_elastic(2,:)
+ b_veloc_elastic(2,:) = 0._CUSTOM_REAL
+ b_accel_elastic(3,:) = b_accel_elastic(2,:)
+ b_accel_elastic(2,:) = 0._CUSTOM_REAL
else !SH (membrane) waves
- do j=1,npoin
- read(55) b_displ_elastic(2,j), &
- b_veloc_elastic(2,j), &
- b_accel_elastic(2,j)
- enddo
- b_displ_elastic(1,:) = ZERO
- b_displ_elastic(3,:) = ZERO
- b_veloc_elastic(1,:) = ZERO
- b_veloc_elastic(3,:) = ZERO
- b_accel_elastic(1,:) = ZERO
- b_accel_elastic(3,:) = ZERO
+ do j=1,npoin
+ read(55) b_displ_elastic(2,j), &
+ b_veloc_elastic(2,j), &
+ b_accel_elastic(2,j)
+ enddo
+ b_displ_elastic(1,:) = 0._CUSTOM_REAL
+ b_displ_elastic(3,:) = 0._CUSTOM_REAL
+ b_veloc_elastic(1,:) = 0._CUSTOM_REAL
+ b_veloc_elastic(3,:) = 0._CUSTOM_REAL
+ b_accel_elastic(1,:) = 0._CUSTOM_REAL
+ b_accel_elastic(3,:) = 0._CUSTOM_REAL
endif
- close(55)
+ close(55)
endif
! poroelastic medium
@@ -5932,20 +5889,23 @@
!---- compute kinetic and potential energy
if(OUTPUT_ENERGY) &
- call compute_energy(displ_elastic,veloc_elastic,displs_poroelastic,velocs_poroelastic, &
+ call 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,poroelastcoef,density, &
+ nspec,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ assign_external_model,it,deltat,t0,kmato,poroelastcoef,density, &
porosity,tortuosity, &
vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext, &
anisotropic,anisotropy,wxgll,wzgll,numat, &
pressure_element,vector_field_element,e1,e11, &
potential_dot_acoustic,potential_dot_dot_acoustic, &
- TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS,p_sv)
!---- display time step and max of norm of displacement
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5 .or. it == NSTEP) then
- call check_stability(myrank,time,it,NSTEP,npoin, &
+ call check_stability(myrank,time,it,NSTEP, &
+ npoin_acoustic,npoin_elastic,npoin_poroelastic, &
any_elastic_glob,any_elastic,displ_elastic, &
any_poroelastic_glob,any_poroelastic, &
displs_poroelastic,displw_poroelastic, &
@@ -5965,7 +5925,8 @@
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,&
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec, &
+ npoin_acoustic,npoin_elastic,npoin_poroelastic,assign_external_model, &
numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext, &
c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,ispec,e1,e11, &
TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
@@ -5974,22 +5935,33 @@
! for acoustic medium, compute vector field from gradient of potential for seismograms
if(seismotype == 1) then
- call compute_vector_one_element(vector_field_element,potential_acoustic,displ_elastic,displs_poroelastic,&
- elastic,poroelastic,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,&
- density,rhoext,assign_external_model)
+ call compute_vector_one_element(vector_field_element,potential_acoustic, &
+ displ_elastic,displs_poroelastic,&
+ elastic,poroelastic,xix,xiz,gammax,gammaz, &
+ ibool,hprime_xx,hprime_zz, &
+ nspec,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ ispec,numat,kmato,density,rhoext,assign_external_model)
else if(seismotype == 2) then
- call compute_vector_one_element(vector_field_element,potential_dot_acoustic,veloc_elastic,velocs_poroelastic, &
- elastic,poroelastic,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,&
- density,rhoext,assign_external_model)
+ call compute_vector_one_element(vector_field_element,potential_dot_acoustic, &
+ veloc_elastic,velocs_poroelastic, &
+ elastic,poroelastic,xix,xiz,gammax,gammaz, &
+ ibool,hprime_xx,hprime_zz, &
+ nspec,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ ispec,numat,kmato,density,rhoext,assign_external_model)
else if(seismotype == 3) then
- call compute_vector_one_element(vector_field_element,potential_dot_dot_acoustic,accel_elastic,accels_poroelastic, &
- elastic,poroelastic,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,&
- density,rhoext,assign_external_model)
+ call compute_vector_one_element(vector_field_element,potential_dot_dot_acoustic, &
+ accel_elastic,accels_poroelastic, &
+ elastic,poroelastic,xix,xiz,gammax,gammaz, &
+ ibool,hprime_xx,hprime_zz, &
+ nspec,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ ispec,numat,kmato,density,rhoext,assign_external_model)
endif
else if(seismotype == 5) then
- call compute_curl_one_element(curl_element,displ_elastic,displs_poroelastic,elastic,poroelastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin, ispec)
+ call compute_curl_one_element(curl_element,displ_elastic, &
+ displs_poroelastic,elastic,poroelastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin_elastic,npoin_poroelastic,ispec)
endif
! perform the general interpolation using Lagrange polynomials
@@ -6441,8 +6413,10 @@
if (myrank == 0) write(IOUT,*) 'drawing displacement vector as small arrows...'
call compute_vector_whole_medium(potential_acoustic,displ_elastic,displs_poroelastic,&
- elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
+ elastic,poroelastic,vector_field_display, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ numat,kmato,density,rhoext,assign_external_model)
call plotpost(vector_field_display,coord,vpext,x_source,z_source,x_final_receiver,z_final_receiver, &
it,deltat,coorg,xinterp,zinterp,shape2D_display, &
@@ -6460,16 +6434,22 @@
myrank,nproc,ier,&
d1_coorg_send_ps_velocity_model,d2_coorg_send_ps_velocity_model, &
d1_coorg_recv_ps_velocity_model,d2_coorg_recv_ps_velocity_model, &
- d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model,d1_RGB_recv_ps_velocity_model,d2_RGB_recv_ps_velocity_model, &
- coorg_send_ps_velocity_model,RGB_send_ps_velocity_model,coorg_recv_ps_velocity_model,RGB_recv_ps_velocity_model, &
- d1_coorg_send_ps_element_mesh,d2_coorg_send_ps_element_mesh,d1_coorg_recv_ps_element_mesh,d2_coorg_recv_ps_element_mesh, &
+ d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model, &
+ d1_RGB_recv_ps_velocity_model,d2_RGB_recv_ps_velocity_model, &
+ coorg_send_ps_velocity_model,RGB_send_ps_velocity_model, &
+ coorg_recv_ps_velocity_model,RGB_recv_ps_velocity_model, &
+ d1_coorg_send_ps_element_mesh,d2_coorg_send_ps_element_mesh, &
+ d1_coorg_recv_ps_element_mesh,d2_coorg_recv_ps_element_mesh, &
d1_color_send_ps_element_mesh,d1_color_recv_ps_element_mesh, &
- coorg_send_ps_element_mesh,color_send_ps_element_mesh,coorg_recv_ps_element_mesh,color_recv_ps_element_mesh, &
+ coorg_send_ps_element_mesh,color_send_ps_element_mesh, &
+ coorg_recv_ps_element_mesh,color_recv_ps_element_mesh, &
d1_coorg_send_ps_abs,d1_coorg_recv_ps_abs,d2_coorg_send_ps_abs,d2_coorg_recv_ps_abs, &
coorg_send_ps_abs,coorg_recv_ps_abs, &
- d1_coorg_send_ps_free_surface,d1_coorg_recv_ps_free_surface,d2_coorg_send_ps_free_surface,d2_coorg_recv_ps_free_surface, &
+ d1_coorg_send_ps_free_surface,d1_coorg_recv_ps_free_surface, &
+ d2_coorg_send_ps_free_surface,d2_coorg_recv_ps_free_surface, &
coorg_send_ps_free_surface,coorg_recv_ps_free_surface, &
- d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field,d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
+ d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field, &
+ d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
coorg_send_ps_vector_field,coorg_recv_ps_vector_field)
else if(imagetype == 2 .and. p_sv) then
@@ -6477,8 +6457,10 @@
if (myrank == 0) write(IOUT,*) 'drawing velocity vector as small arrows...'
call compute_vector_whole_medium(potential_dot_acoustic,veloc_elastic,velocs_poroelastic,&
- elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
+ elastic,poroelastic,vector_field_display, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ numat,kmato,density,rhoext,assign_external_model)
call plotpost(vector_field_display,coord,vpext,x_source,z_source,x_final_receiver,z_final_receiver, &
it,deltat,coorg,xinterp,zinterp,shape2D_display, &
@@ -6496,16 +6478,22 @@
myrank,nproc,ier,&
d1_coorg_send_ps_velocity_model,d2_coorg_send_ps_velocity_model, &
d1_coorg_recv_ps_velocity_model,d2_coorg_recv_ps_velocity_model, &
- d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model,d1_RGB_recv_ps_velocity_model,d2_RGB_recv_ps_velocity_model, &
- coorg_send_ps_velocity_model,RGB_send_ps_velocity_model,coorg_recv_ps_velocity_model,RGB_recv_ps_velocity_model, &
- d1_coorg_send_ps_element_mesh,d2_coorg_send_ps_element_mesh,d1_coorg_recv_ps_element_mesh,d2_coorg_recv_ps_element_mesh, &
+ d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model, &
+ d1_RGB_recv_ps_velocity_model,d2_RGB_recv_ps_velocity_model, &
+ coorg_send_ps_velocity_model,RGB_send_ps_velocity_model, &
+ coorg_recv_ps_velocity_model,RGB_recv_ps_velocity_model, &
+ d1_coorg_send_ps_element_mesh,d2_coorg_send_ps_element_mesh, &
+ d1_coorg_recv_ps_element_mesh,d2_coorg_recv_ps_element_mesh, &
d1_color_send_ps_element_mesh,d1_color_recv_ps_element_mesh, &
- coorg_send_ps_element_mesh,color_send_ps_element_mesh,coorg_recv_ps_element_mesh,color_recv_ps_element_mesh, &
+ coorg_send_ps_element_mesh,color_send_ps_element_mesh, &
+ coorg_recv_ps_element_mesh,color_recv_ps_element_mesh, &
d1_coorg_send_ps_abs,d1_coorg_recv_ps_abs,d2_coorg_send_ps_abs,d2_coorg_recv_ps_abs, &
coorg_send_ps_abs,coorg_recv_ps_abs, &
- d1_coorg_send_ps_free_surface,d1_coorg_recv_ps_free_surface,d2_coorg_send_ps_free_surface,d2_coorg_recv_ps_free_surface, &
+ d1_coorg_send_ps_free_surface,d1_coorg_recv_ps_free_surface, &
+ d2_coorg_send_ps_free_surface,d2_coorg_recv_ps_free_surface, &
coorg_send_ps_free_surface,coorg_recv_ps_free_surface, &
- d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field,d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
+ d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field, &
+ d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
coorg_send_ps_vector_field,coorg_recv_ps_vector_field)
else if(imagetype == 3 .and. p_sv) then
@@ -6513,8 +6501,10 @@
if (myrank == 0) write(IOUT,*) 'drawing acceleration vector as small arrows...'
call compute_vector_whole_medium(potential_dot_dot_acoustic,accel_elastic,accels_poroelastic,&
- elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
+ elastic,poroelastic,vector_field_display, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ numat,kmato,density,rhoext,assign_external_model)
call plotpost(vector_field_display,coord,vpext,x_source,z_source,x_final_receiver,z_final_receiver, &
it,deltat,coorg,xinterp,zinterp,shape2D_display, &
@@ -6532,16 +6522,22 @@
myrank,nproc,ier,&
d1_coorg_send_ps_velocity_model,d2_coorg_send_ps_velocity_model, &
d1_coorg_recv_ps_velocity_model,d2_coorg_recv_ps_velocity_model, &
- d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model,d1_RGB_recv_ps_velocity_model,d2_RGB_recv_ps_velocity_model, &
- coorg_send_ps_velocity_model,RGB_send_ps_velocity_model,coorg_recv_ps_velocity_model,RGB_recv_ps_velocity_model, &
- d1_coorg_send_ps_element_mesh,d2_coorg_send_ps_element_mesh,d1_coorg_recv_ps_element_mesh,d2_coorg_recv_ps_element_mesh, &
+ d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model, &
+ d1_RGB_recv_ps_velocity_model,d2_RGB_recv_ps_velocity_model, &
+ coorg_send_ps_velocity_model,RGB_send_ps_velocity_model, &
+ coorg_recv_ps_velocity_model,RGB_recv_ps_velocity_model, &
+ d1_coorg_send_ps_element_mesh,d2_coorg_send_ps_element_mesh, &
+ d1_coorg_recv_ps_element_mesh,d2_coorg_recv_ps_element_mesh, &
d1_color_send_ps_element_mesh,d1_color_recv_ps_element_mesh, &
- coorg_send_ps_element_mesh,color_send_ps_element_mesh,coorg_recv_ps_element_mesh,color_recv_ps_element_mesh, &
+ coorg_send_ps_element_mesh,color_send_ps_element_mesh, &
+ coorg_recv_ps_element_mesh,color_recv_ps_element_mesh, &
d1_coorg_send_ps_abs,d1_coorg_recv_ps_abs,d2_coorg_send_ps_abs,d2_coorg_recv_ps_abs, &
coorg_send_ps_abs,coorg_recv_ps_abs, &
- d1_coorg_send_ps_free_surface,d1_coorg_recv_ps_free_surface,d2_coorg_send_ps_free_surface,d2_coorg_recv_ps_free_surface, &
+ d1_coorg_send_ps_free_surface,d1_coorg_recv_ps_free_surface, &
+ d2_coorg_send_ps_free_surface,d2_coorg_recv_ps_free_surface, &
coorg_send_ps_free_surface,coorg_recv_ps_free_surface, &
- d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field,d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
+ d1_coorg_send_ps_vector_field,d1_coorg_recv_ps_vector_field, &
+ d2_coorg_send_ps_vector_field,d2_coorg_recv_ps_vector_field, &
coorg_send_ps_vector_field,coorg_recv_ps_vector_field)
else if(imagetype == 4 .or. .not. p_sv) then
@@ -6568,24 +6564,30 @@
if (myrank == 0) write(IOUT,*) 'drawing image of z (if P-SV) or y (if SH) component of displacement vector...'
call compute_vector_whole_medium(potential_acoustic,displ_elastic,displs_poroelastic,&
- elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
+ elastic,poroelastic,vector_field_display, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ numat,kmato,density,rhoext,assign_external_model)
else if(imagetype == 2) then
if (myrank == 0) write(IOUT,*) 'drawing image of z (if P-SV) or y (if SH) component of velocity vector...'
call compute_vector_whole_medium(potential_dot_acoustic,veloc_elastic,velocs_poroelastic,&
- elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
+ elastic,poroelastic,vector_field_display, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ numat,kmato,density,rhoext,assign_external_model)
else if(imagetype == 3) then
if (myrank == 0) write(IOUT,*) 'drawing image of z (if P-SV) or y (if SH) component of acceleration vector...'
call compute_vector_whole_medium(potential_dot_dot_acoustic,accel_elastic,accels_poroelastic,&
- elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
+ elastic,poroelastic,vector_field_display, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz, &
+ nspec,npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic, &
+ numat,kmato,density,rhoext,assign_external_model)
else if(imagetype == 4 .and. p_sv) then
@@ -6593,7 +6595,8 @@
call compute_pressure_whole_medium(potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec, &
+ npoin,npoin_acoustic,npoin_elastic,npoin_poroelastic,assign_external_model, &
numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext, &
c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,e1,e11, &
TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
@@ -6706,26 +6709,26 @@
!--- save last frame
!
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1 .and. any_elastic) then
- if ( myrank == 0 ) then
- write(IOUT,*)
- write(IOUT,*) 'Saving elastic last frame...'
- write(IOUT,*)
- endif
+ if ( myrank == 0 ) then
+ write(IOUT,*)
+ write(IOUT,*) 'Saving elastic last frame...'
+ write(IOUT,*)
+ endif
write(outputname,'(a,i6.6,a)') 'lastframe_elastic',myrank,'.bin'
open(unit=55,file='OUTPUT_FILES/'//outputname,status='unknown',form='unformatted')
- if(p_sv)then !P-SV waves
- do j=1,npoin
- write(55) displ_elastic(1,j), displ_elastic(3,j), &
+ if(p_sv)then !P-SV waves
+ do j=1,npoin
+ write(55) displ_elastic(1,j), displ_elastic(3,j), &
veloc_elastic(1,j), veloc_elastic(3,j), &
accel_elastic(1,j), accel_elastic(3,j)
- enddo
- else !SH (membrane) waves
- do j=1,npoin
- write(55) displ_elastic(2,j), &
+ enddo
+ else !SH (membrane) waves
+ do j=1,npoin
+ write(55) displ_elastic(2,j), &
veloc_elastic(2,j), &
accel_elastic(2,j)
- enddo
- endif
+ enddo
+ endif
close(55)
endif
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2011-02-23 01:47:34 UTC (rev 17948)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2011-02-23 03:14:03 UTC (rev 17949)
@@ -341,7 +341,7 @@
write(11,110) FACTORXSU,NSTEP,deltat,-t0,nrec
do irec=1,nrec
- !daniel
+ ! this format statement might now work for larger meshes
!write(11,"(f12.5)") st_xval(irec)
write(11,*) st_xval(irec)
if(irec < nrec) write(11,*) ','
More information about the CIG-COMMITS
mailing list