[cig-commits] r12970 - seismo/3D/SPECFEM3D/branches/update_temporary
nlegoff at geodynamics.org
nlegoff at geodynamics.org
Thu Sep 25 06:58:04 PDT 2008
Author: nlegoff
Date: 2008-09-25 06:58:03 -0700 (Thu, 25 Sep 2008)
New Revision: 12970
Modified:
seismo/3D/SPECFEM3D/branches/update_temporary/assemble_MPI_scalar.f90
seismo/3D/SPECFEM3D/branches/update_temporary/constants.h.in
seismo/3D/SPECFEM3D/branches/update_temporary/locate_receivers.f90
seismo/3D/SPECFEM3D/branches/update_temporary/meshfem3D.f90
seismo/3D/SPECFEM3D/branches/update_temporary/parallel.f90
seismo/3D/SPECFEM3D/branches/update_temporary/specfem3D.f90
Log:
added locating receivers on the surface and recording seismos normal to surface in case of external mesh.
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/assemble_MPI_scalar.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/assemble_MPI_scalar.f90 2008-09-25 03:43:14 UTC (rev 12969)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/assemble_MPI_scalar.f90 2008-09-25 13:58:03 UTC (rev 12970)
@@ -277,3 +277,85 @@
endif
end subroutine assemble_MPI_scalar_ext_mesh
+
+!
+!----
+!
+
+ subroutine assemble_MPI_scalar_i_ext_mesh(NPROC,NGLOB_AB,array_val, &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh &
+ )
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+! array to assemble
+ integer, dimension(NGLOB_AB) :: array_val
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+
+ integer, dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh
+
+ integer :: ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh
+
+ integer ipoin,iinterface
+ integer sender,receiver
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! partition border copy into the buffer
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ buffer_send_scalar_ext_mesh(ipoin,iinterface) = array_val(ibool_interfaces_ext_mesh(ipoin,iinterface))
+ enddo
+ enddo
+
+! send messages
+ do iinterface = 1, ninterfaces_ext_mesh
+ call issend_i(buffer_send_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_scalar_ext_mesh(iinterface) &
+ )
+ call irecv_i(buffer_recv_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_scalar_ext_mesh(iinterface) &
+ )
+ enddo
+
+! wait for communications completion
+ do iinterface = 1, ninterfaces_ext_mesh
+ call wait_req(request_recv_scalar_ext_mesh(iinterface))
+ enddo
+
+! adding contributions of neighbours
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_scalar_ext_mesh(ipoin,iinterface)
+ enddo
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_scalar_i_ext_mesh
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/constants.h.in 2008-09-25 03:43:14 UTC (rev 12969)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/constants.h.in 2008-09-25 13:58:03 UTC (rev 12970)
@@ -122,6 +122,12 @@
! whether or not an external mesh is used (provided by CUBIT for example)
logical, parameter :: USE_EXTERNAL_MESH = .false.
+! no lagrange interpolation on seismograms (we take the value on one NGLL point)
+ logical, parameter :: FASTER_RECEIVERS_POINTS_ONLY = .false.
+
+! the receivers can be located inside the model
+ logical, parameter :: RECEIVERS_CAN_BE_BURIED_EXT_MESH = .true.
+
! deltat
double precision, parameter :: DT_ext_mesh = 0.001d0
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/locate_receivers.f90 2008-09-25 03:43:14 UTC (rev 12969)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/locate_receivers.f90 2008-09-25 13:58:03 UTC (rev 12970)
@@ -32,7 +32,9 @@
xi_receiver,eta_receiver,gamma_receiver,station_name,network_name,nu, &
NPROC,utm_x_source,utm_y_source, &
TOPOGRAPHY,itopo_bathy,UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
- NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO)
+ NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
+ iglob_is_surface_external_mesh,ispec_is_surface_external_mesh &
+ )
implicit none
@@ -57,6 +59,12 @@
integer itopo_bathy(NX_TOPO,NY_TOPO)
double precision long_corner,lat_corner,ratio_xi,ratio_eta
+! for surface locating and normal computing with external mesh
+ integer :: pt0_ix,pt0_iy,pt0_iz,pt1_ix,pt1_iy,pt1_iz,pt2_ix,pt2_iy,pt2_iz
+ real(kind=CUSTOM_REAL), dimension(3) :: u_vector,v_vector,w_vector
+ logical, dimension(NGLOB_AB) :: iglob_is_surface_external_mesh
+ logical, dimension(NSPEC_AB) :: ispec_is_surface_external_mesh
+
integer, allocatable, dimension(:) :: ix_initial_guess,iy_initial_guess,iz_initial_guess
integer iprocloop
@@ -68,7 +76,7 @@
double precision, allocatable, dimension(:,:) :: x_found_all,y_found_all,z_found_all
integer irec
- integer i,j,k,ispec,iglob
+ integer i,j,k,ispec,iglob,imin,imax,jmin,jmax,kmin,kmax
integer icornerlong,icornerlat
double precision utm_x_source,utm_y_source
@@ -120,7 +128,8 @@
integer, allocatable, dimension(:,:) :: ispec_selected_rec_all
double precision, allocatable, dimension(:) :: stlat,stlon,stele,stbur,stutm_x,stutm_y
double precision, allocatable, dimension(:,:) :: xi_receiver_all,eta_receiver_all,gamma_receiver_all
-
+ double precision, allocatable, dimension(:,:,:,:) :: nu_all
+
character(len=150) OUTPUT_FILES
! **************
@@ -179,6 +188,7 @@
allocate(y_found_all(nrec,0:NPROC-1))
allocate(z_found_all(nrec,0:NPROC-1))
allocate(final_distance_all(nrec,0:NPROC-1))
+ allocate(nu_all(3,3,nrec,0:NPROC-1))
! loop on all the stations
do irec=1,nrec
@@ -190,7 +200,8 @@
if (ios /= 0) call exit_mpi(myrank, 'Error reading station file '//trim(rec_filename))
! convert station location to UTM
- call utm_geo(stlon(irec),stlat(irec),stutm_x(irec),stutm_y(irec),UTM_PROJECTION_ZONE,ILONGLAT2UTM,SUPPRESS_UTM_PROJECTION)
+ call utm_geo(stlon(irec),stlat(irec),stutm_x(irec),stutm_y(irec),UTM_PROJECTION_ZONE,ILONGLAT2UTM, &
+ (SUPPRESS_UTM_PROJECTION .or. USE_EXTERNAL_MESH))
! compute horizontal distance between source and receiver in km
horiz_dist(irec) = dsqrt((stutm_y(irec)-utm_y_source)**2 + (stutm_x(irec)-utm_x_source)**2) / 1000.
@@ -220,6 +231,7 @@
nu(3,2,irec) = 0.d0
nu(3,3,irec) = 1.d0
+ if (.not. USE_EXTERNAL_MESH) then
! compute elevation of topography at the receiver location
! we assume that receivers are always at the surface i.e. not buried
if(TOPOGRAPHY) then
@@ -265,21 +277,56 @@
z_target(irec) = elevation(irec) - stbur(irec)
if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
+ else
+
+ x_target(irec) = stutm_x(irec)
+ y_target(irec) = stutm_y(irec)
+ z_target(irec) = stbur(irec)
+ if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
+
+ endif ! of if (.not. USE_EXTERNAL_MESH)
+
! examine top of the elements only (receivers always at the surface)
! k = NGLLZ
do ispec=1,NSPEC_AB
-! modification by Qinya Liu: idoubling is no longer used because receivers can now be in depth
-! (in previous versions, receivers were always assumed to be at the surface)
+! define the interval in which we look for points
+ if(FASTER_RECEIVERS_POINTS_ONLY) then
+ imin = 1
+ imax = NGLLX
+ jmin = 1
+ jmax = NGLLY
+
+ kmin = 1
+ kmax = NGLLZ
+
+ else
! loop only on points inside the element
! exclude edges to ensure this point is not shared with other elements
- do k=2,NGLLZ-1
- do j=2,NGLLY-1
- do i=2,NGLLX-1
+ imin = 2
+ imax = NGLLX - 1
+ jmin = 2
+ jmax = NGLLY - 1
+
+ kmin = 2
+ kmax = NGLLZ - 1
+ endif
+
+ do k = kmin,kmax
+ do j = jmin,jmax
+ do i = imin,imax
+
iglob = ibool(i,j,k,ispec)
+
+ if (USE_EXTERNAL_MESH .and. (.not. RECEIVERS_CAN_BE_BURIED_EXT_MESH)) then
+ if ((.not. iglob_is_surface_external_mesh(iglob)) .or. (.not. ispec_is_surface_external_mesh(ispec))) then
+ cycle
+ endif
+ endif
+
dist = dsqrt((x_target(irec)-dble(xstore(iglob)))**2 &
+(y_target(irec)-dble(ystore(iglob)))**2 &
+(z_target(irec)-dble(zstore(iglob)))**2)
@@ -291,16 +338,163 @@
ix_initial_guess(irec) = i
iy_initial_guess(irec) = j
iz_initial_guess(irec) = k
+
+ xi_receiver(irec) = dble(ix_initial_guess(irec))
+ eta_receiver(irec) = dble(iy_initial_guess(irec))
+ gamma_receiver(irec) = dble(iz_initial_guess(irec))
+ x_found(irec) = xstore(iglob)
+ y_found(irec) = ystore(iglob)
+ z_found(irec) = zstore(iglob)
endif
enddo
enddo
enddo
+
+! compute final distance between asked and found (converted to km)
+ final_distance(irec) = dsqrt((x_target(irec)-x_found(irec))**2 + &
+ (y_target(irec)-y_found(irec))**2 + (z_target(irec)-z_found(irec))**2)
! endif
! end of loop on all the spectral elements in current slice
enddo
+! get normal to the face of the hexaedra if receiver is on the surface
+ if (USE_EXTERNAL_MESH .and. (.not. RECEIVERS_CAN_BE_BURIED_EXT_MESH)) then
+ pt0_ix = -1
+ pt0_iy = -1
+ pt0_iz = -1
+ pt1_ix = -1
+ pt1_iy = -1
+ pt1_iz = -1
+ pt2_ix = -1
+ pt2_iy = -1
+ pt2_iz = -1
+! we get two vectors of the face (three points) to compute the normal
+ if (ix_initial_guess(irec) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(1,2,2,ispec_selected_rec(irec)))) then
+ pt0_ix = 1
+ pt0_iy = NGLLY
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = 1
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+ if (ix_initial_guess(irec) == NGLLX .and. &
+ iglob_is_surface_external_mesh(ibool(NGLLX,2,2,ispec_selected_rec(irec)))) then
+ pt0_ix = NGLLX
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = NGLLX
+ pt1_iy = NGLLY
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = 1
+ pt2_iz = NGLLZ
+ endif
+ if (iy_initial_guess(irec) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(2,1,2,ispec_selected_rec(irec)))) then
+ pt0_ix = 1
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = NGLLX
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = 1
+ pt2_iy = 1
+ pt2_iz = NGLLZ
+ endif
+ if (iy_initial_guess(irec) == NGLLY .and. &
+ iglob_is_surface_external_mesh(ibool(2,NGLLY,2,ispec_selected_rec(irec)))) then
+ pt0_ix = NGLLX
+ pt0_iy = NGLLY
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = NGLLY
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+ if (iz_initial_guess(irec) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(2,2,1,ispec_selected_rec(irec)))) then
+ pt0_ix = NGLLX
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = NGLLY
+ pt2_iz = 1
+ endif
+ if (iz_initial_guess(irec) == NGLLZ .and. &
+ iglob_is_surface_external_mesh(ibool(2,2,NGLLZ,ispec_selected_rec(irec)))) then
+ pt0_ix = 1
+ pt0_iy = 1
+ pt0_iz = NGLLZ
+ pt1_ix = NGLLX
+ pt1_iy = 1
+ pt1_iz = NGLLZ
+ pt2_ix = 1
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+
+ if (pt0_ix<0 .or.pt0_iy<0 .or. pt0_iz<0 .or. &
+ pt1_ix<0 .or. pt1_iy<0 .or. pt1_iz<0 .or. &
+ pt2_ix<0 .or. pt2_iy<0 .or. pt2_iz<0) then
+ stop 'error in computing normal for receivers.'
+ endif
+
+ u_vector(1) = xstore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_rec(irec))) &
+ - xstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ u_vector(2) = ystore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_rec(irec))) &
+ - ystore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ u_vector(3) = zstore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_rec(irec))) &
+ - zstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ v_vector(1) = xstore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_rec(irec))) &
+ - xstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ v_vector(2) = ystore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_rec(irec))) &
+ - ystore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ v_vector(3) = zstore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_rec(irec))) &
+ - zstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+
+! cross product
+ w_vector(1) = u_vector(2)*v_vector(3) - u_vector(3)*v_vector(2)
+ w_vector(2) = u_vector(3)*v_vector(1) - u_vector(1)*v_vector(3)
+ w_vector(3) = u_vector(1)*v_vector(2) - u_vector(2)*v_vector(1)
+
+! normalize vector w
+ w_vector(:) = w_vector(:)/sqrt(w_vector(1)**2+w_vector(2)**2+w_vector(3)**2)
+
+! build the two other vectors for a direct base : we normalize u, and v=w^u
+ u_vector(:) = u_vector(:)/sqrt(u_vector(1)**2+u_vector(2)**2+u_vector(3)**2)
+ v_vector(1) = w_vector(2)*u_vector(3) - w_vector(3)*u_vector(2)
+ v_vector(2) = w_vector(3)*u_vector(1) - w_vector(1)*u_vector(3)
+ v_vector(3) = w_vector(1)*u_vector(2) - w_vector(2)*u_vector(1)
+
+! build rotation matrice nu for seismograms
+! East (u)
+ nu(1,1,irec) = u_vector(1)
+ nu(1,2,irec) = v_vector(1)
+ nu(1,3,irec) = w_vector(1)
+
+! North (v)
+ nu(2,1,irec) = u_vector(2)
+ nu(2,2,irec) = v_vector(2)
+ nu(2,3,irec) = w_vector(2)
+
+! Vertical (w)
+ nu(3,1,irec) = u_vector(3)
+ nu(3,2,irec) = v_vector(3)
+ nu(3,3,irec) = w_vector(3)
+
+ endif ! of if (USE_EXTERNAL_MESH .and. (.not. RECEIVERS_CAN_BE_BURIED_EXT_MESH))
+
! end of loop on all the stations
enddo
@@ -311,6 +505,8 @@
! find the best (xi,eta,gamma) for each receiver
! ****************************************
+ if(.not. FASTER_RECEIVERS_POINTS_ONLY) then
+
! loop on all the receivers to iterate in that slice
do irec = 1,nrec
@@ -424,6 +620,8 @@
enddo
+ endif ! of if (.not. FASTER_RECEIVERS_POINTS_ONLY)
+
! synchronize all the processes to make sure all the estimates are available
call sync_all()
@@ -437,6 +635,7 @@
call gather_all_dp(x_found,nrec,x_found_all,nrec,NPROC)
call gather_all_dp(y_found,nrec,y_found_all,nrec,NPROC)
call gather_all_dp(z_found,nrec,z_found_all,nrec,NPROC)
+ call gather_all_dp(nu,3*3*nrec,nu_all,3*3*nrec,NPROC)
! this is executed by main process only
if(myrank == 0) then
@@ -459,6 +658,7 @@
x_found(irec) = x_found_all(irec,iprocloop)
y_found(irec) = y_found_all(irec,iprocloop)
z_found(irec) = z_found_all(irec,iprocloop)
+ nu(:,:,irec) = nu_all(:,:,irec,iprocloop)
endif
enddo
final_distance(irec) = distmin
@@ -481,7 +681,15 @@
write(IMAIN,*) 'closest estimate found: ',sngl(final_distance(irec)),' m away'
write(IMAIN,*) ' in slice ',islice_selected_rec(irec),' in element ',ispec_selected_rec(irec)
- write(IMAIN,*) ' at xi,eta,gamma coordinates = ',xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec)
+ if(FASTER_RECEIVERS_POINTS_ONLY) then
+ write(IMAIN,*) 'in point i,j,k = ',nint(xi_receiver(irec)),nint(eta_receiver(irec)),nint(gamma_receiver(irec))
+ !write(IMAIN,*) 'in point i,j,k = ',x_found(irec),y_found(irec),z_found(irec)
+ write(IMAIN,*) 'nu1 = ',nu(1,:,irec)
+ write(IMAIN,*) 'nu2 = ',nu(2,:,irec)
+ write(IMAIN,*) 'nu3 = ',nu(3,:,irec)
+ else
+ write(IMAIN,*) ' at xi,eta,gamma coordinates = ',xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec)
+ endif
! add warning if estimate is poor
! (usually means receiver outside the mesh given by the user)
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/meshfem3D.f90 2008-09-25 03:43:14 UTC (rev 12969)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/meshfem3D.f90 2008-09-25 13:58:03 UTC (rev 12970)
@@ -920,7 +920,8 @@
open(unit=IIN,file=rec_filename,status='old',action='read')
do irec = 1,nrec
read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
- if(stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ if((stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ .or. USE_EXTERNAL_MESH) &
nrec_filtered = nrec_filtered + 1
enddo
close(IIN)
@@ -938,7 +939,8 @@
do irec = 1,nrec
read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
- if(stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ if((stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ .or. USE_EXTERNAL_MESH) &
write(IOUT,*) station_name(1:len_trim(station_name)),' ',network_name(1:len_trim(network_name)),' ', &
sngl(stlat),' ',sngl(stlon), ' ', sngl(stele), ' ', sngl(stbur)
enddo
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/parallel.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/parallel.f90 2008-09-25 03:43:14 UTC (rev 12969)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/parallel.f90 2008-09-25 13:58:03 UTC (rev 12970)
@@ -460,6 +460,60 @@
!----
!
+ subroutine issend_i(sendbuf, sendcount, dest, sendtag, req)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+ integer sendcount, dest, sendtag, req
+ integer, dimension(sendcount) :: sendbuf
+
+! MPI status of messages to be received
+ integer msg_status(MPI_STATUS_SIZE)
+
+ integer ier
+
+ call MPI_ISSEND(sendbuf(1),sendcount,MPI_INTEGER,dest,sendtag, &
+ MPI_COMM_WORLD,req,ier)
+
+ end subroutine issend_i
+
+!
+!----
+!
+
+ subroutine irecv_i(recvbuf, recvcount, dest, recvtag, req)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+ integer recvcount, dest, recvtag, req
+ integer, dimension(recvcount) :: recvbuf
+
+! MPI status of messages to be received
+ integer msg_status(MPI_STATUS_SIZE)
+
+ integer ier
+
+ call MPI_IRECV(recvbuf(1),recvcount,MPI_INTEGER,dest,recvtag, &
+ MPI_COMM_WORLD,req,ier)
+
+ end subroutine irecv_i
+
+!
+!----
+!
+
subroutine wait_req(req)
implicit none
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/specfem3D.f90 2008-09-25 03:43:14 UTC (rev 12969)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/specfem3D.f90 2008-09-25 13:58:03 UTC (rev 12970)
@@ -454,6 +454,14 @@
integer, dimension(:), allocatable :: request_send_vector_ext_mesh
integer, dimension(:), allocatable :: request_recv_vector_ext_mesh
+! for detecting surface receivers and source in case of external mesh
+ integer, dimension(:), allocatable :: valence_external_mesh
+ logical, dimension(:), allocatable :: iglob_is_surface_external_mesh
+ logical, dimension(:), allocatable :: ispec_is_surface_external_mesh
+ integer, dimension(:,:), allocatable :: buffer_send_scalar_i_ext_mesh
+ integer, dimension(:,:), allocatable :: buffer_recv_scalar_i_ext_mesh
+ integer :: ii,jj,kk
+
! ************** PROGRAM STARTS HERE **************
! sizeprocs returns number of processes started
@@ -884,6 +892,82 @@
endif ! end of (.not. USE_EXTERNAL_MESH)
+! detecting surface points/elements (based on valence check on NGLL points) for external mesh
+ allocate(valence_external_mesh(NGLOB_AB))
+ allocate(ispec_is_surface_external_mesh(NSPEC_AB))
+ allocate(iglob_is_surface_external_mesh(NGLOB_AB))
+
+ if (.not. RECEIVERS_CAN_BE_BURIED_EXT_MESH) then
+ valence_external_mesh(:) = 0
+ ispec_is_surface_external_mesh(:) = .false.
+ iglob_is_surface_external_mesh(:) = .false.
+ do ispec = 1, NSPEC_AB
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ valence_external_mesh(iglob) = valence_external_mesh(iglob) + 1
+
+ enddo
+ enddo
+ enddo
+
+ enddo
+
+ allocate(buffer_send_scalar_i_ext_mesh(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+ allocate(buffer_recv_scalar_i_ext_mesh(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+
+ call assemble_MPI_scalar_i_ext_mesh(NPROC,NGLOB_AB,valence_external_mesh, &
+ buffer_send_scalar_i_ext_mesh,buffer_recv_scalar_i_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh &
+ )
+
+ do ispec = 1, NSPEC_AB
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ if ( &
+ (k == 1 .or. k == NGLLZ) .and. (j /= 1 .and. j /= NGLLY) .and. (i /= 1 .and. i /= NGLLX) .or. &
+ (j == 1 .or. j == NGLLY) .and. (k /= 1 .and. k /= NGLLZ) .and. (i /= 1 .and. i /= NGLLX) .or. &
+ (i == 1 .or. i == NGLLX) .and. (k /= 1 .and. k /= NGLLZ) .and. (j /= 1 .and. j /= NGLLY) &
+ ) then
+ iglob = ibool(i,j,k,ispec)
+ if (valence_external_mesh(iglob) == 1) then
+ ispec_is_surface_external_mesh(ispec) = .true.
+ if (k == 1 .or. k == NGLLZ) then
+ do jj = 1, NGLLY
+ do ii = 1, NGLLX
+ iglob_is_surface_external_mesh(ibool(ii,jj,k,ispec)) = .true.
+ enddo
+ enddo
+ endif
+ if (j == 1 .or. j == NGLLY) then
+ do kk = 1, NGLLZ
+ do ii = 1, NGLLX
+ iglob_is_surface_external_mesh(ibool(ii,j,kk,ispec)) = .true.
+ enddo
+ enddo
+ endif
+ if (i == 1 .or. i == NGLLX) then
+ do kk = 1, NGLLZ
+ do jj = 1, NGLLY
+ iglob_is_surface_external_mesh(ibool(i,jj,kk,ispec)) = .true.
+ enddo
+ enddo
+ endif
+ endif
+
+ endif
+ enddo
+ enddo
+ enddo
+
+ enddo
+
+ endif !
+
! $$$$$$$$$$$$$$$$$$$$$$$$ SOURCES $$$$$$$$$$$$$$$$$
! read topography and bathymetry file
@@ -1021,7 +1105,9 @@
xi_receiver,eta_receiver,gamma_receiver,station_name,network_name,nu, &
NPROC,utm_x_source(1),utm_y_source(1), &
TOPOGRAPHY,itopo_bathy,UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
- NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO)
+ NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
+ iglob_is_surface_external_mesh,ispec_is_surface_external_mesh &
+)
!###################### SOURCE ARRAYS ################
@@ -2721,6 +2807,22 @@
irec = number_receiver_global(irec_local)
! perform the general interpolation using Lagrange polynomials
+ if(FASTER_RECEIVERS_POINTS_ONLY) then
+
+ iglob = ibool(nint(xi_receiver(irec)),nint(eta_receiver(irec)), &
+ nint(gamma_receiver(irec)),ispec_selected_rec(irec))
+ dxd = dble(displ(1,iglob))
+ dyd = dble(displ(2,iglob))
+ dzd = dble(displ(3,iglob))
+ vxd = dble(veloc(1,iglob))
+ vyd = dble(veloc(2,iglob))
+ vzd = dble(veloc(3,iglob))
+ axd = dble(accel(1,iglob))
+ ayd = dble(accel(2,iglob))
+ azd = dble(accel(3,iglob))
+
+ else
+
dxd = ZERO
dyd = ZERO
dzd = ZERO
@@ -2763,6 +2865,8 @@
enddo
enddo
+
+
else if (SIMULATION_TYPE == 2) then
do k = 1,NGLLZ
@@ -2832,8 +2936,8 @@
enddo
endif
+ endif ! end of if(FASTER_RECEIVERS_POINTS_ONLY)
-
! store North, East and Vertical components
! distinguish between single and double precision for reals
More information about the cig-commits
mailing list