[cig-commits] r19349 - in seismo/3D/SPECFEM3D/trunk: src/decompose_mesh_SCOTCH src/specfem3D utils/adjoint_sources
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Jan 11 14:32:26 PST 2012
Author: danielpeter
Date: 2012-01-11 14:32:26 -0800 (Wed, 11 Jan 2012)
New Revision: 19349
Added:
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/SU_adjoint.f90
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
Log:
adds source encoding capabilities for acoustic sources; updates receiver location handling for SU_FORMAT outputs
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2012-01-11 05:52:51 UTC (rev 19348)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2012-01-11 22:32:26 UTC (rev 19349)
@@ -1424,6 +1424,10 @@
endif
enddo
+ ! checks if any poroelastic/elastic elements are set
+ !if( .not. any(is_poroelastic) ) return
+ !if( .not. any(is_elastic) ) return
+
! gets neighbors by 4 common nodes (face)
allocate(xadj(0:nelmnts),stat=ier)
if( ier /= 0 ) stop 'error allocating array xadj'
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2012-01-11 05:52:51 UTC (rev 19348)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2012-01-11 22:32:26 UTC (rev 19349)
@@ -39,7 +39,8 @@
use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
- station_name,network_name,adj_source_file,nrec_local,number_receiver_global
+ station_name,network_name,adj_source_file,nrec_local,number_receiver_global, &
+ pm1_source_encoding
implicit none
include "constants.h"
@@ -145,6 +146,9 @@
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ ! source encoding
+ stf_used = stf_used * pm1_source_encoding(isource)
+
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
! to add minus the source to Chi_dot_dot to get plus the source in pressure:
@@ -164,6 +168,9 @@
! quasi-heaviside
!stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ ! source encoding
+ stf = stf * pm1_source_encoding(isource)
+
! distinguishes between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
stf_used = sngl(stf)
@@ -382,6 +389,9 @@
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
+ ! source encoding
+ stf_used = stf_used * pm1_source_encoding(isource)
+
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
! to add minus the source to Chi_dot_dot to get plus the source in pressure:
@@ -398,6 +408,12 @@
! gaussian source time
stf = comp_source_time_function_gauss(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ ! quasi-heaviside
+ !stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+
+ ! source encoding
+ stf = stf * pm1_source_encoding(isource)
+
! distinguishes between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
stf_used = sngl(stf)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2012-01-11 05:52:51 UTC (rev 19348)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2012-01-11 22:32:26 UTC (rev 19349)
@@ -40,30 +40,48 @@
include "constants.h"
- logical SUPPRESS_UTM_PROJECTION
+ integer :: myrank
+ integer :: NSPEC_AB,NGLOB_AB
- integer NPROC,UTM_PROJECTION_ZONE
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
- integer nrec,myrank
+ ! arrays containing coordinates of the points
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
- integer NSPEC_AB,NGLOB_AB
+ ! Gauss-Lobatto-Legendre points of integration
+ double precision xigll(NGLLX)
+ double precision yigll(NGLLY)
+ double precision zigll(NGLLZ)
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ ! input receiver file name
+ character(len=*) rec_filename
-! arrays containing coordinates of the points
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
+ ! receivers
+ integer :: nrec
+ integer, dimension(nrec),intent(out) :: islice_selected_rec,ispec_selected_rec
+ double precision, dimension(nrec),intent(out) :: xi_receiver,eta_receiver,gamma_receiver
+ double precision, dimension(3,3,nrec),intent(out) :: nu
+ character(len=MAX_LENGTH_STATION_NAME), dimension(nrec),intent(out) :: station_name
+ character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec),intent(out) :: network_name
-! 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
+ integer :: NPROC,UTM_PROJECTION_ZONE
+ double precision :: utm_x_source,utm_y_source
+ logical :: SUPPRESS_UTM_PROJECTION
+
+ ! free surface
integer :: num_free_surface_faces
- 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, dimension(num_free_surface_faces) :: free_surface_ispec
integer, dimension(3,NGLLSQUARE,num_free_surface_faces) :: free_surface_ijk
+ ! local parameters
+
+ ! for surface locating and normal computing with external mesh
+ real(kind=CUSTOM_REAL), dimension(3) :: u_vector,v_vector,w_vector
+ integer :: pt0_ix,pt0_iy,pt0_iz,pt1_ix,pt1_iy,pt1_iz,pt2_ix,pt2_iy,pt2_iz
+
integer, allocatable, dimension(:) :: ix_initial_guess,iy_initial_guess,iz_initial_guess
-
integer iprocloop
integer ios
@@ -73,64 +91,45 @@
double precision, allocatable, dimension(:) :: x_target,y_target,z_target
double precision, allocatable, dimension(:) :: horiz_dist
double precision, allocatable, dimension(:) :: x_found,y_found,z_found
+ double precision dist
+ double precision xi,eta,gamma,dx,dy,dz,dxi,deta,dgamma
+ double precision x,y,z
+ double precision xix,xiy,xiz
+ double precision etax,etay,etaz
+ double precision gammax,gammay,gammaz
+ ! coordinates of the control points of the surface element
+ double precision xelm(NGNOD),yelm(NGNOD),zelm(NGNOD)
+
integer irec
integer i,j,k,ispec,iglob,iface,inode,imin,imax,jmin,jmax,kmin,kmax,igll,jgll,kgll
integer iselected,jselected,iface_selected,iadjust,jadjust
integer iproc(1)
- double precision utm_x_source,utm_y_source
- double precision dist
- double precision xi,eta,gamma,dx,dy,dz,dxi,deta,dgamma
-
-! Gauss-Lobatto-Legendre points of integration
- double precision xigll(NGLLX)
- double precision yigll(NGLLY)
- double precision zigll(NGLLZ)
-
-! input receiver file name
- character(len=*) rec_filename
-
-! topology of the control points of the surface element
+ ! topology of the control points of the surface element
integer iax,iay,iaz
integer iaddx(NGNOD),iaddy(NGNOD),iaddz(NGNOD)
-! coordinates of the control points of the surface element
- double precision xelm(NGNOD),yelm(NGNOD),zelm(NGNOD)
-
integer iter_loop,ispec_iterate
-
integer ia
- double precision x,y,z
- double precision xix,xiy,xiz
- double precision etax,etay,etaz
- double precision gammax,gammay,gammaz
-! timer MPI
+ ! timer MPI
double precision, external :: wtime
double precision time_start,tCPU
-! use dynamic allocation
+ ! use dynamic allocation
double precision, dimension(:), allocatable :: final_distance
double precision distmin,final_distance_max
-! receiver information
-! timing information for the stations
-! station information for writing the seismograms
-
+ ! receiver information
+ ! station information for writing the seismograms
integer :: iglob_selected
- integer, dimension(nrec) :: islice_selected_rec,ispec_selected_rec
- double precision, dimension(nrec) :: xi_receiver,eta_receiver,gamma_receiver
- double precision, dimension(3,3,nrec) :: nu
- character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
- character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
double precision, allocatable, dimension(:) :: stlat,stlon,stele,stbur,stutm_x,stutm_y,elevation
-
double precision, allocatable, dimension(:) :: x_found_all,y_found_all,z_found_all
double precision, dimension(:), allocatable :: final_distance_all
- integer, allocatable, dimension(:) :: ispec_selected_rec_all
double precision, allocatable, dimension(:) :: xi_receiver_all,eta_receiver_all,gamma_receiver_all
double precision, allocatable, dimension(:,:,:) :: nu_all
+ integer, allocatable, dimension(:) :: ispec_selected_rec_all
integer :: ier
character(len=256) OUTPUT_FILES
@@ -140,12 +139,30 @@
integer :: imin_temp,imax_temp,jmin_temp,jmax_temp,kmin_temp,kmax_temp
integer,dimension(NGLLX*NGLLY*NGLLZ) :: iglob_temp
-! **************
+ ! SU_FORMAT parameters
+ double precision :: llat,llon,lele,lbur
+ logical :: SU_station_file_exists
+ ! get MPI starting time
+ time_start = wtime()
+
+ ! user output
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) '********************'
+ write(IMAIN,*) ' locating receivers'
+ write(IMAIN,*) '********************'
+ write(IMAIN,*)
+ write(IMAIN,'(1x,a,a,a)') 'reading receiver information from ', trim(rec_filename), ' file'
+ write(IMAIN,*)
+ endif
+
! dimension of model in current proc
xmin=minval(xstore(:)); xmax=maxval(xstore(:))
ymin=minval(ystore(:)); ymax=maxval(ystore(:))
zmin=minval(zstore(:)); zmax=maxval(zstore(:))
+
+ ! loop limits
if(FASTER_RECEIVERS_POINTS_ONLY) then
imin_temp = 1; imax_temp = NGLLX
jmin_temp = 1; jmax_temp = NGLLY
@@ -156,31 +173,73 @@
kmin_temp = 2; kmax_temp = NGLLZ - 1
endif
- ! get MPI starting time
- time_start = wtime()
-
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) '********************'
- write(IMAIN,*) ' locating receivers'
- write(IMAIN,*) '********************'
- write(IMAIN,*)
- endif
-
! define topology of the control element
call usual_hex_nodes(iaddx,iaddy,iaddz)
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) '*****************************************************************'
- write(IMAIN,'(1x,a,a,a)') 'reading receiver information from ', trim(rec_filename), ' file'
- write(IMAIN,*) '*****************************************************************'
- endif
-
- ! get number of stations from receiver file
+ ! opens STATIONS file
open(unit=1,file=trim(rec_filename),status='old',action='read',iostat=ios)
if (ios /= 0) call exit_mpi(myrank,'error opening file '//trim(rec_filename))
+ ! get the base pathname for output files
+ call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', OUTPUT_FILES_PATH(1:len_trim(OUTPUT_FILES_PATH)))
+
+ ! checks if station locations already available
+ if( SU_FORMAT ) then
+ ! checks if file with station infos located from previous run exists
+ INQUIRE(file=trim(OUTPUT_FILES)//'/SU_stations_info.bin',exist=SU_station_file_exists)
+ if ( SU_station_file_exists ) then
+ ! all processes read in stations names from STATIONS file
+ do irec=1,nrec
+ read(1,*,iostat=ios) station_name(irec),network_name(irec),llat,llon,lele,lbur
+ if (ios /= 0) call exit_mpi(myrank, 'Error reading station file '//trim(rec_filename))
+ enddo
+ close(1)
+ ! master reads in available station information
+ if( myrank == 0 ) then
+ open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/SU_stations_info.bin', &
+ status='old',action='read',form='unformatted',iostat=ios)
+ if (ios /= 0) call exit_mpi(myrank,'error opening file '//trim(rec_filename))
+ write(IMAIN,*) 'station details from SU_stations_info.bin'
+ allocate(x_found(nrec),y_found(nrec),z_found(nrec))
+ ! reads in station infos
+ read(IOUT_SU) islice_selected_rec,ispec_selected_rec
+ read(IOUT_SU) xi_receiver,eta_receiver,gamma_receiver
+ read(IOUT_SU) x_found,y_found,z_found
+ read(IOUT_SU) nu
+ close(IOUT_SU)
+ ! write the locations of stations, so that we can load them and write them to SU headers later
+ open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/output_list_stations.txt', &
+ status='unknown',action='write',iostat=ios)
+ if( ios /= 0 ) call exit_mpi(myrank,'error opening file '//trim(OUTPUT_FILES)//'/output_list_stations.txt')
+ do irec=1,nrec
+ write(IOUT_SU,*) x_found(irec),y_found(irec),z_found(irec)
+ enddo
+ close(IOUT_SU)
+ deallocate(x_found,y_found,z_found)
+ endif
+ ! main process broadcasts the results to all the slices
+ call bcast_all_i(islice_selected_rec,nrec)
+ call bcast_all_i(ispec_selected_rec,nrec)
+ call bcast_all_dp(xi_receiver,nrec)
+ call bcast_all_dp(eta_receiver,nrec)
+ call bcast_all_dp(gamma_receiver,nrec)
+ call bcast_all_dp(nu,NDIM*NDIM*nrec)
+ call sync_all()
+ ! user output
+ if( myrank == 0 ) then
+ ! elapsed time since beginning of mesh generation
+ tCPU = wtime() - time_start
+ write(IMAIN,*)
+ write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',tCPU
+ write(IMAIN,*)
+ write(IMAIN,*) 'End of receiver detection - done'
+ write(IMAIN,*)
+ endif
+ ! everything done
+ return
+ endif
+ endif
+
! allocate memory for arrays using number of stations
allocate(stlat(nrec), &
stlon(nrec), &
@@ -225,11 +284,15 @@
horiz_dist(irec) = dsqrt((stutm_y(irec)-utm_y_source)**2 + (stutm_x(irec)-utm_x_source)**2) / 1000.
! print some information about stations
- if(myrank == 0) &
+ if(myrank == 0) then
+ ! limits user output if too many receivers
+ if( nrec < 1000 .and. ( .not. SU_FORMAT ) ) then
write(IMAIN,*) 'Station #',irec,': ',station_name(irec)(1:len_trim(station_name(irec))), &
'.',network_name(irec)(1:len_trim(network_name(irec))), &
' horizontal distance: ',sngl(horiz_dist(irec)),' km'
-
+ endif
+ endif
+
! get approximate topography elevation at source long/lat coordinates
! set distance to huge initial value
distmin = HUGEVAL
@@ -255,8 +318,8 @@
iglob = ibool(igll,jgll,kgll,ispec)
! keep this point if it is closer to the receiver
- dist = dsqrt((stutm_x(irec)-dble(xstore(iglob)))**2 + &
- (stutm_y(irec)-dble(ystore(iglob)))**2)
+ dist = (stutm_x(irec)-dble(xstore(iglob)))**2 + &
+ (stutm_y(irec)-dble(ystore(iglob)))**2
if(dist < distmin) then
distmin = dist
iglob_selected = iglob
@@ -285,7 +348,7 @@
elevation_node(inode) = zstore(iglob)
dist_node(inode) = dsqrt((stutm_x(irec)-dble(xstore(iglob)))**2 + &
- (stutm_y(irec)-dble(ystore(iglob)))**2)
+ (stutm_y(irec)-dble(ystore(iglob)))**2)
inode = inode + 1
end do
end do
@@ -293,9 +356,9 @@
if(dist < distmin) then
distmin = dist
altitude_rec(1) = (dist_node(1)/dist)*elevation_node(1) + &
- (dist_node(2)/dist)*elevation_node(2) + &
- (dist_node(3)/dist)*elevation_node(3) + &
- (dist_node(4)/dist)*elevation_node(4)
+ (dist_node(2)/dist)*elevation_node(2) + &
+ (dist_node(3)/dist)*elevation_node(3) + &
+ (dist_node(4)/dist)*elevation_node(4)
endif
end do
end do
@@ -344,127 +407,125 @@
!if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
if (.not. SU_FORMAT) then
- ! determines closest GLL point
- ispec_selected_rec(irec) = 0
- do ispec=1,NSPEC_AB
+ ! determines closest GLL point
+ ispec_selected_rec(irec) = 0
+ do ispec=1,NSPEC_AB
- ! define the interval in which we look for points
- if(FASTER_RECEIVERS_POINTS_ONLY) then
- imin = 1
- imax = NGLLX
+ ! 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
+ imin = 2
+ imax = NGLLX - 1
+ jmin = 2
+ jmax = NGLLY - 1
+ kmin = 2
+ kmax = NGLLZ - 1
+ endif
- jmin = 1
- jmax = NGLLY
+ do k = kmin,kmax
+ do j = jmin,jmax
+ do i = imin,imax
- kmin = 1
- kmax = NGLLZ
+ iglob = ibool(i,j,k,ispec)
- else
- ! loop only on points inside the element
- ! exclude edges to ensure this point is not shared with other elements
- imin = 2
- imax = NGLLX - 1
+ if (.not. RECVS_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
- jmin = 2
- jmax = NGLLY - 1
+ dist = (x_target(irec)-dble(xstore(iglob)))**2 &
+ + (y_target(irec)-dble(ystore(iglob)))**2 &
+ + (z_target(irec)-dble(zstore(iglob)))**2
- kmin = 2
- kmax = NGLLZ - 1
- 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(irec) = i
+ iy_initial_guess(irec) = j
+ iz_initial_guess(irec) = k
- do k = kmin,kmax
- do j = jmin,jmax
- do i = imin,imax
+ 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
- iglob = ibool(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
- if (.not. RECVS_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
+ ! 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)
- dist = dsqrt((x_target(irec)-dble(xstore(iglob)))**2 &
- +(y_target(irec)-dble(ystore(iglob)))**2 &
- +(z_target(irec)-dble(zstore(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(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)
-
- ! end of loop on all the spectral elements in current slice
- enddo
+ ! end of loop on all the spectral elements in current slice
+ enddo
else
- ispec_selected_rec(irec) = 0
- ix_initial_guess(irec) = 0
- iy_initial_guess(irec) = 0
- iz_initial_guess(irec) = 0
- final_distance(irec) = HUGEVAL
- if ( (x_target(irec)>=xmin .and. x_target(irec)<=xmax) .and. &
- (y_target(irec)>=ymin .and. y_target(irec)<=ymax) .and. &
- (z_target(irec)>=zmin .and. z_target(irec)<=zmax) ) then
- do ispec=1,NSPEC_AB
- iglob_temp=reshape(ibool(:,:,:,ispec),(/NGLLX*NGLLY*NGLLZ/))
- xmin_ELE=minval(xstore(iglob_temp))
- xmax_ELE=maxval(xstore(iglob_temp))
- ymin_ELE=minval(ystore(iglob_temp))
- ymax_ELE=maxval(ystore(iglob_temp))
- zmin_ELE=minval(zstore(iglob_temp))
- zmax_ELE=maxval(zstore(iglob_temp))
- if ( (x_target(irec)>=xmin_ELE .and. x_target(irec)<=xmax_ELE) .and. &
- (y_target(irec)>=ymin_ELE .and. y_target(irec)<=ymax_ELE) .and. &
- (z_target(irec)>=zmin_ELE .and. z_target(irec)<=zmax_ELE) ) then
- ! we find the element (ispec) which "may" contain the receiver (irec)
- ! so we only need to compute distances (which is expensive because of "dsqrt") within those elements
- ispec_selected_rec(irec) = ispec
- do k = kmin_temp,kmax_temp
- do j = jmin_temp,jmax_temp
- do i = imin_temp,imax_temp
- iglob = ibool(i,j,k,ispec)
- ! for comparison purpose, we don't have to do "dsqrt", which is expensive
- dist = ((x_target(irec)-dble(xstore(iglob)))**2 &
- +(y_target(irec)-dble(ystore(iglob)))**2 &
- +(z_target(irec)-dble(zstore(iglob)))**2)
- if(dist < distmin) then
- distmin = dist
- 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
+ ispec_selected_rec(irec) = 0
+ ix_initial_guess(irec) = 0
+ iy_initial_guess(irec) = 0
+ iz_initial_guess(irec) = 0
+ final_distance(irec) = HUGEVAL
+ if ( (x_target(irec)>=xmin .and. x_target(irec)<=xmax) .and. &
+ (y_target(irec)>=ymin .and. y_target(irec)<=ymax) .and. &
+ (z_target(irec)>=zmin .and. z_target(irec)<=zmax) ) then
+ do ispec=1,NSPEC_AB
+ iglob_temp=reshape(ibool(:,:,:,ispec),(/NGLLX*NGLLY*NGLLZ/))
+ xmin_ELE=minval(xstore(iglob_temp))
+ xmax_ELE=maxval(xstore(iglob_temp))
+ ymin_ELE=minval(ystore(iglob_temp))
+ ymax_ELE=maxval(ystore(iglob_temp))
+ zmin_ELE=minval(zstore(iglob_temp))
+ zmax_ELE=maxval(zstore(iglob_temp))
+ if ( (x_target(irec)>=xmin_ELE .and. x_target(irec)<=xmax_ELE) .and. &
+ (y_target(irec)>=ymin_ELE .and. y_target(irec)<=ymax_ELE) .and. &
+ (z_target(irec)>=zmin_ELE .and. z_target(irec)<=zmax_ELE) ) then
+ ! we find the element (ispec) which "may" contain the receiver (irec)
+ ! so we only need to compute distances (which is expensive because of "dsqrt") within those elements
+ ispec_selected_rec(irec) = ispec
+ do k = kmin_temp,kmax_temp
+ do j = jmin_temp,jmax_temp
+ do i = imin_temp,imax_temp
+ iglob = ibool(i,j,k,ispec)
+ ! for comparison purpose, we don't have to do "dsqrt", which is expensive
+ dist = ((x_target(irec)-dble(xstore(iglob)))**2 &
+ + (y_target(irec)-dble(ystore(iglob)))**2 &
+ + (z_target(irec)-dble(zstore(iglob)))**2)
+ if(dist < distmin) then
+ distmin = dist
+ 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
- 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 ! if receiver "may" be within this element
- enddo ! do ispec=1,NSPEC_AB
- endif ! if receiver "may" be within this proc
+ enddo
+ enddo
+ 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 ! if receiver "may" be within this element
+ enddo ! do ispec=1,NSPEC_AB
+ endif ! if receiver "may" be within this proc
endif !if (.not. SU_FORMAT)
if (ispec_selected_rec(irec) == 0) then
@@ -806,11 +867,18 @@
do irec=1,nrec
+ ! checks stations location
+ if(final_distance(irec) == HUGEVAL) then
+ write(IMAIN,*) 'error locating station # ',irec,' ',station_name(irec),network_name(irec)
+ call exit_MPI(myrank,'error locating receiver')
+ endif
+
+ ! limits user output if too many receivers
+ if( nrec < 1000 .and. (.not. SU_FORMAT ) ) then
+
write(IMAIN,*)
write(IMAIN,*) 'station # ',irec,' ',station_name(irec),network_name(irec)
- if(final_distance(irec) == HUGEVAL) call exit_MPI(myrank,'error locating receiver')
-
write(IMAIN,*) ' original latitude: ',sngl(stlat(irec))
write(IMAIN,*) ' original longitude: ',sngl(stlon(irec))
if( SUPPRESS_UTM_PROJECTION ) then
@@ -856,7 +924,6 @@
endif
write(IMAIN,*)
-
! add warning if estimate is poor
! (usually means receiver outside the mesh given by the user)
if(final_distance(irec) > 3000.d0) then
@@ -867,6 +934,8 @@
write(IMAIN,*)
+ endif
+
enddo
! compute maximal distance for all the receivers
@@ -886,9 +955,6 @@
write(IMAIN,*) '************************************************************'
endif
- ! get the base pathname for output files
- call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', OUTPUT_FILES_PATH(1:len_trim(OUTPUT_FILES_PATH)))
-
!! write the list of stations and associated epicentral distance
!open(unit=27,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
!do irec=1,nrec
@@ -897,12 +963,27 @@
!close(27)
! write the locations of stations, so that we can load them and write them to SU headers later
- open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
+ open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/output_list_stations.txt', &
+ status='unknown',action='write',iostat=ios)
+ if( ios /= 0 ) call exit_mpi(myrank,'error opening file '//trim(OUTPUT_FILES)//'/output_list_stations.txt')
do irec=1,nrec
write(IOUT_SU,*) x_found(irec),y_found(irec),z_found(irec)
enddo
close(IOUT_SU)
+ ! stores station infos for later runs
+ if( SU_FORMAT ) then
+ open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/SU_stations_info.bin', &
+ status='unknown',action='write',form='unformatted',iostat=ios)
+ if( ios == 0 ) then
+ write(IOUT_SU) islice_selected_rec,ispec_selected_rec
+ write(IOUT_SU) xi_receiver,eta_receiver,gamma_receiver
+ write(IOUT_SU) x_found,y_found,z_found
+ write(IOUT_SU) nu
+ close(IOUT_SU)
+ endif
+ endif
+
! elapsed time since beginning of mesh generation
tCPU = wtime() - time_start
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-01-11 05:52:51 UTC (rev 19348)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-01-11 22:32:26 UTC (rev 19349)
@@ -95,6 +95,10 @@
nu_source(3,3,NSOURCES),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for sources'
+ ! for source encoding (acoustic sources so far only)
+ allocate(pm1_source_encoding(NSOURCES),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for sources'
+
! locate sources in the mesh
!
! returns: islice_selected_source & ispec_selected_source,
@@ -487,12 +491,15 @@
! user output
call any_all_l( is_on, is_on_all )
if( myrank == 0 .and. is_on_all ) then
- write(IMAIN,*) '**********************************************************************'
- write(IMAIN,*) '*** station:',irec,' ***'
- write(IMAIN,*) '*** Warning: acoustic receiver located exactly on the free surface ***'
- write(IMAIN,*) '*** Warning: tangential component will be zero there ***'
- write(IMAIN,*) '**********************************************************************'
- write(IMAIN,*)
+ ! limits user output if too many receivers
+ if( nrec < 1000 .and. (.not. SU_FORMAT ) ) then
+ write(IMAIN,*) '**********************************************************************'
+ write(IMAIN,*) '*** station:',irec,' ***'
+ write(IMAIN,*) '*** Warning: acoustic receiver located exactly on the free surface ***'
+ write(IMAIN,*) '*** Warning: tangential component will be zero there ***'
+ write(IMAIN,*) '**********************************************************************'
+ write(IMAIN,*)
+ endif
endif
enddo ! num_free_surface_faces
@@ -559,6 +566,12 @@
! (and getting rid of 1/sqrt(2) factor from scalar moment tensor definition above)
factor_source = factor_source * sqrt(2.0) / sqrt(3.0)
+ ! source encoding
+ ! determines factor +/-1 depending on sign of moment tensor
+ ! (see e.g. Krebs et al., 2009. Fast full-wavefield seismic inversion using encoded sources,
+ ! Geophysics, 74 (6), WCC177-WCC188.)
+ pm1_source_encoding(isource) = sign(1.0,Mxx(isource))
+
! source array interpolated on all element gll points (only used for non point sources)
call compute_arrays_source_acoustic(xi_source(isource),eta_source(isource),gamma_source(isource),&
sourcearray,xigll,yigll,zigll,factor_source)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-01-11 05:52:51 UTC (rev 19348)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-01-11 22:32:26 UTC (rev 19349)
@@ -98,6 +98,10 @@
double precision :: t0
real(kind=CUSTOM_REAL) :: stf_used_total
integer :: NSOURCES
+ ! source encoding
+ ! for acoustic sources: takes +/- 1 sign, depending on sign(Mxx)[ = sign(Myy) = sign(Mzz)
+ ! since they have to equal in the acoustic setting]
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: pm1_source_encoding
! receiver information
character(len=256) :: rec_filename,filtered_rec_filename,dummystring
@@ -355,7 +359,7 @@
module specfem_par_acoustic
-! parameter module for elastic solver
+! parameter module for acoustic solver
use constants,only: CUSTOM_REAL
implicit none
Added: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/SU_adjoint.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/SU_adjoint.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/SU_adjoint.f90 2012-01-11 22:32:26 UTC (rev 19349)
@@ -0,0 +1,197 @@
+!
+! program to create adjoint source files in SU (Seismic Unix) format
+!
+! this template program uses SU-formatted files.
+! by default, it takes a pure waveform misfit to construct adjoint sources (adj = syn - dat)
+! please modify routine create_adjoint_trace() according to your needs.
+!
+! for example, compile with:
+! > gfortran -o ../../example/homogeneous/halfspace/bin/xSU_adjoint SU_adjoint.f90
+!
+! call in example directory example/homogeneous_halfspace/ by:
+! > ./bin/xSU_adjoint 4
+!
+!-------------------------------------------------------------------------------------------------
+
+ module SU_adjoint_var
+
+ implicit none
+
+ real(kind=4),dimension(:),allocatable :: syn,dat,adj
+
+ character(len=128) :: DATA_PATH = "./DATA/"
+ character(len=128) :: SYN_PATH = "./in_out_files/OUTPUT_FILES/"
+ character(len=128) :: ADJ_PATH = "./in_out_files/SEM/"
+
+ end module
+
+!-------------------------------------------------------------------------------------------------
+
+ program SU_adjoint
+
+ use SU_adjoint_var
+ implicit none
+
+ integer :: NPROC
+ integer :: iproc,icomp,ios,irec,i
+
+ character(len=512) :: filename
+ character(len=512) :: arg(1)
+ character(len=16) :: procname
+ character(len=1),dimension(3) :: compstr = (/ "x" , "y" , "z" /)
+
+ ! SU header
+ integer(kind=4) :: r4head(60)
+ integer(kind=4) :: header4(1)
+ integer(kind=2) :: header2(2)
+ equivalence(header2,header4)
+
+ integer :: NSTEP
+ double precision :: DT
+
+ ! reads in file arguments
+ i = 1
+ do while (1 == 1)
+ call getarg(i,arg(i))
+ if (i <= 1 .and. trim(arg(i)) == '') then
+ print*,'Usage: '
+ print*,' ./xSU_adjoint NPROC'
+ print*,'with'
+ print*,' NPROC: total number of partitions'
+ stop
+ endif
+ if (trim(arg(i)) == '') exit
+ if (i == 1) then
+ read(arg(i),*,iostat=ios) NPROC
+ if (ios /= 0) stop 'Error reading NPROC'
+ endif
+ i = i+1
+ enddo
+
+ ! user output
+ print*, "SU adjoint "
+ print*
+ print*, "number of partitions : ",NPROC
+ print*, "data path : ",trim(DATA_PATH)
+ print*, "synthetics path : ",trim(SYN_PATH)
+ print*, "adjoint sources path : ",trim(ADJ_PATH)
+ print*
+
+ ! reads NSTEP info from seismograms
+ ! only do this once
+ iproc = 0
+ write(procname,"(i4)") iproc
+ filename = trim(DATA_PATH)//trim(adjustl(procname))//"_dx_SU"
+ open(11,file=trim(filename),access='direct',status='old', &
+ recl=240,iostat=ios)
+ if( ios /= 0 ) then
+ print*,'error opening file: ',trim(filename)
+ stop 'error opening data file'
+ endif
+
+ ! reads in header
+ read(11,rec=1,iostat=ios) r4head
+ close(11)
+
+ ! extracts header information
+ header4 = r4head(29)
+ NSTEP = header2(2)
+ header4 = r4head(30)
+ DT = header2(1) * 1.0d-6
+
+ ! user output
+ print*, "NSTEP = ",NSTEP
+ print*, "DT = ",DT
+ print*
+
+ ! allocates arrays
+ allocate(syn(NSTEP),dat(NSTEP),adj(NSTEP))
+
+ ! loops over all partitions
+ do iproc = 0, NPROC - 1
+
+ ! partition name
+ write(procname,"(i4)") iproc
+ print*,"...reading partition: ",iproc
+
+ ! loops over components
+ do icomp = 1,3
+ ! user output
+ print*," component ",compstr(icomp)
+
+ ! opens files
+ ! data
+ filename = trim(DATA_PATH)//trim(adjustl(procname))//"_d"//compstr(icomp)//"_SU"
+ open(11,file=trim(filename),access='direct',status='old', &
+ recl=240+4*NSTEP,iostat=ios)
+ if( ios /= 0 ) then
+ print*,'error opening file: ',trim(filename)
+ stop 'error opening input data file '
+ endif
+
+ ! synthetics
+ filename = trim(SYN_PATH)//trim(adjustl(procname))//"_d"//compstr(icomp)//"_SU"
+ open(22,file=trim(filename),access='direct',status='old', &
+ recl=240+4*NSTEP,iostat=ios)
+ if( ios /= 0 ) then
+ print*,'error opening file: ',trim(filename)
+ stop 'error opening input file '
+ endif
+
+ ! adjoint sources
+ filename = trim(ADJ_PATH)//trim(adjustl(procname))//"_d"//compstr(icomp)//"_SU"//".adj"
+ open(33,file=trim(filename),access='direct',status='unknown', &
+ recl=240+4*NSTEP,iostat = ios)
+ if( ios /= 0 ) then
+ print*,'error opening file: ',trim(filename)
+ stop 'error opening output file '
+ endif
+
+ ! loops over all records
+ irec=1
+ do while(ios==0)
+ ! reads in data
+ read(11,rec=irec,iostat=ios) r4head,dat
+ if (ios /= 0) cycle
+
+ ! reads in synthetics
+ read(22,rec=irec,iostat=ios) r4head,syn
+ if (ios /= 0) cycle
+
+ ! creates adjoint trace
+ call create_adjoint_trace()
+
+ ! writes out adjoint source
+ write(33,rec=irec,iostat=ios) r4head,adj
+ if (ios /= 0) cycle
+
+ irec=irec+1
+ enddo
+ close(11)
+ close(22)
+ close(33)
+ enddo
+
+ ! user output
+ print*, " receivers read: ",irec
+ enddo
+
+ ! user output
+ print*, "done adjoint sources"
+ print*, "please, check output files in directory: ",trim(ADJ_PATH)
+
+ end program
+
+!-------------------------------------------------------------------------------------------------
+
+ subroutine create_adjoint_trace()
+
+ use SU_adjoint_var
+ implicit none
+
+ ! MODIFY THIS according to your misfit function
+
+ ! waveform misfit:
+ adj(:) = syn(:) - dat(:)
+
+ end subroutine
More information about the CIG-COMMITS
mailing list