[cig-commits] r19248 - in seismo/3D/SPECFEM3D/trunk: doc/USER_MANUAL src/shared src/specfem3D utils
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Nov 29 07:09:45 PST 2011
Author: danielpeter
Date: 2011-11-29 07:09:45 -0800 (Tue, 29 Nov 2011)
New Revision: 19248
Added:
seismo/3D/SPECFEM3D/trunk/utils/locate_partition.f90
Modified:
seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.pdf
seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex
seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
Log:
adds normalization option to create_movie_shakemap_AVS_DX_GMT.f90; updates movie outputs for saving displacements; updates manual in movies section
Modified: seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex
===================================================================
--- seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex 2011-11-29 04:46:30 UTC (rev 19247)
+++ seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex 2011-11-29 15:09:45 UTC (rev 19248)
@@ -1747,7 +1747,8 @@
\texttt{MOVIE\_VOLUME}, produces large output files. \texttt{MOVIE\_VOLUME}
files are saved in the \texttt{LOCAL\_PATH} directory, whereas \texttt{MOVIE\_SURFACE}
output files are saved in the \texttt{in\_out\_files/OUTPUT\_FILES} directory. We
-save the velocity field. The look of a movie is determined by the
+save the displacement field if the parameter \texttt{SAVE\_DISPLACEMENT} is set, otherwise the velocity field is saved.
+The look of a movie is determined by the
half-duration of the source. The half-duration should be large enough
so that the movie does not contain frequencies that are not resolved
by the mesh, i.e., it should not contain numerical noise. This can
@@ -1761,10 +1762,10 @@
\sqrt{(}\mathrm{\mathtt{HALF\_DURATIO}\mathtt{N}^{2}}+\mathrm{\mathtt{HDUR\_MOVI}\mathtt{E}^{2}})\]
\textbf{NOTE:} If \texttt{HDUR\_MOVIE} is set to 0.0, the code will
select the appropriate value of 1.1 $\times$ smallest period. As
-usual, for a point source one can set \texttt{HALF\_DURATION} in the
-\texttt{Par\_file} to be 0.0 and \texttt{HDUR\_MOVIE} = 0.0 to get
+usual, for a point source one can set \texttt{half duration} in the
+\texttt{CMTSOLUTION} file to be 0.0 and \texttt{HDUR\_MOVIE} = 0.0 in the \texttt{Par\_file} to get
the highest frequencies resolved by the simulation, but for a finite
-source one would keep all the \texttt{HALF\_DURATIONs} as prescribed
+source one would keep all the \texttt{half durations} as prescribed
by the finite source model and set \texttt{HDUR\_MOVIE} = 0.0.
\end{quote}
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2011-11-29 04:46:30 UTC (rev 19247)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2011-11-29 15:09:45 UTC (rev 19248)
@@ -35,10 +35,14 @@
implicit none
include "constants.h"
- include "surface_from_mesher.h"
+ include "../../in_out_files/OUTPUT_FILES/surface_from_mesher.h"
!-------------------------------------------------------------------------------------------------
! user parameters
+
+! normalizes field display values
+ logical, parameter :: NORMALIZE_OUTPUT = .false.
+
! threshold in percent of the maximum below which we cut the amplitude
logical, parameter :: APPLY_THRESHOLD = .false.
real(kind=CUSTOM_REAL), parameter :: THRESHOLD = 1._CUSTOM_REAL / 100._CUSTOM_REAL
@@ -548,44 +552,44 @@
min_field_current = minval(field_display(:))
max_field_current = maxval(field_display(:))
- ! print minimum and maximum amplitude in current snapshot
- print *
- print *,'minimum amplitude in current snapshot = ',min_field_current
- print *,'maximum amplitude in current snapshot = ',max_field_current
- print *
-
if(plot_shaking_map) then
- ! compute min and max of data value to normalize
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
! print minimum and maximum amplitude in current snapshot
print *
print *,'minimum amplitude in current snapshot after removal = ',min_field_current
print *,'maximum amplitude in current snapshot after removal = ',max_field_current
print *
+ else
+ ! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot = ',min_field_current
+ print *,'maximum amplitude in current snapshot = ',max_field_current
+ print *
endif
! apply scaling in all cases for movies
if(.not. plot_shaking_map) then
- ! make sure range is always symmetric and center is in zero
- ! this assumption works only for fields that can be negative
- ! would not work for norm of vector for instance
- ! (we would lose half of the color palette if no negative values)
- max_absol = max(abs(min_field_current),abs(max_field_current))
- min_field_current = - max_absol
- max_field_current = + max_absol
+ ! normalizes values
+ if( NORMALIZE_OUTPUT ) then
+ ! make sure range is always symmetric and center is in zero
+ ! this assumption works only for fields that can be negative
+ ! would not work for norm of vector for instance
+ ! (we would lose half of the color palette if no negative values)
+ max_absol = max(abs(min_field_current),abs(max_field_current))
+ min_field_current = - max_absol
+ max_field_current = + max_absol
- ! normalize field to [0:1]
- if( abs(max_field_current - min_field_current) > TINYVAL ) &
- field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
+ ! normalize field to [0:1]
+ if( abs(max_field_current - min_field_current) > TINYVAL ) &
+ field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
- ! rescale to [-1,1]
- field_display(:) = 2.*field_display(:) - 1.
+ ! rescale to [-1,1]
+ field_display(:) = 2.*field_display(:) - 1.
- ! apply threshold to normalized field
- if(APPLY_THRESHOLD) &
- where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
+ ! apply threshold to normalized field
+ if(APPLY_THRESHOLD) &
+ where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
+ endif
! apply non linear scaling to normalized field if needed
if(NONLINEAR_SCALING) then
@@ -596,13 +600,15 @@
endwhere
endif
- ! map back to [0,1]
- field_display(:) = (field_display(:) + 1.) / 2.
+ ! normalizes values
+ if( NORMALIZE_OUTPUT ) then
+ ! map back to [0,1]
+ field_display(:) = (field_display(:) + 1.) / 2.
- ! map field to [0:255] for AVS color scale
- field_display(:) = 255. * field_display(:)
+ ! map field to [0:255] for AVS color scale
+ field_display(:) = 255. * field_display(:)
+ endif
-
! apply scaling only if selected for shaking map
else if(NONLINEAR_SCALING .and. iscaling_shake == 1) then
@@ -668,21 +674,21 @@
if(USE_GMT) then
-! output list of points
- mask_point = .false.
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
- ibool_number = iglob(ilocnum+ieoff)
- if(.not. mask_point(ibool_number)) then
- call utm_geo(long,lat,xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff), &
- UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
- write(11,*) long,lat,field_display(ilocnum+ieoff)
- endif
- mask_point(ibool_number) = .true.
- enddo
- enddo
+ ! output list of points
+ mask_point = .false.
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ do ilocnum = 1,NGNOD2D_AVS_DX
+ ibool_number = iglob(ilocnum+ieoff)
+ if(.not. mask_point(ibool_number)) then
+ call utm_geo(long,lat,xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff), &
+ UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
+ write(11,*) long,lat,field_display(ilocnum+ieoff)
+ endif
+ mask_point(ibool_number) = .true.
+ enddo
+ enddo
else
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2011-11-29 04:46:30 UTC (rev 19247)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2011-11-29 15:09:45 UTC (rev 19248)
@@ -329,12 +329,12 @@
use specfem_par_movie
implicit none
- real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: veloc_element
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: val_element
real(kind=CUSTOM_REAL),dimension(1):: dummy
integer :: ispec2D,ispec,ipoin,iglob,ier
! allocate array for single elements
- allocate( veloc_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ allocate( val_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for movie elements'
! initializes arrays for point coordinates
@@ -365,51 +365,40 @@
ispec = faces_surface_ext_mesh_ispec(ispec2D)
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
- potential_dot_acoustic, veloc_element,&
+ if(SAVE_DISPLACEMENT) then
+ ! displacement vector
+ call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
+ potential_acoustic, val_element,&
hprime_xx,hprime_yy,hprime_zz, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
ibool,rhostore)
+ else
+ ! velocity vector
+ call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
+ potential_dot_acoustic, val_element,&
+ hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ ibool,rhostore)
+ endif
endif
if (USE_HIGHRES_FOR_MOVIES) then
+ ! all GLL points on a face
do ipoin = 1, NGLLX*NGLLY
iglob = faces_surface_ext_mesh(ipoin,ispec2D)
- ! saves velocity vector
- if( ispec_is_elastic(ispec) ) then
- ! velocity x,y,z-components
- store_val_ux_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = veloc(1,iglob)
- store_val_uy_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = veloc(2,iglob)
- store_val_uz_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = veloc(3,iglob)
- endif
-
- ! acoustic pressure potential
- if( ispec_is_acoustic(ispec) ) then
- ! puts velocity values into storage array
- call wmo_get_vel_vector(ispec,ispec2D,ipoin, &
- veloc_element, &
+ ! puts displ/velocity values into storage array
+ call wmo_get_vel_vector(ispec,ispec2D,ipoin,iglob, &
+ val_element, &
NGLLX*NGLLY)
- endif
enddo
else
+ ! only corner points
do ipoin = 1, 4
iglob = faces_surface_ext_mesh(ipoin,ispec2D)
- ! saves velocity vector
- if( ispec_is_elastic(ispec) ) then
- ! velocity x,y,z-components
- store_val_ux_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = veloc(1,iglob)
- store_val_uy_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = veloc(2,iglob)
- store_val_uz_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = veloc(3,iglob)
- endif
-
- ! acoustic pressure potential
- if( ispec_is_acoustic(ispec) ) then
- ! puts velocity values into storage array
- call wmo_get_vel_vector(ispec,ispec2D,ipoin, &
- veloc_element, &
+ ! puts displ/velocity values into storage array
+ call wmo_get_vel_vector(ispec,ispec2D,ipoin,iglob, &
+ val_element, &
NGNOD2D)
- endif
enddo
endif
enddo
@@ -480,45 +469,65 @@
close(IOUT)
endif
- deallocate(veloc_element)
+ deallocate(val_element)
end subroutine wmo_create_movie_surface_em
!================================================================
- subroutine wmo_get_vel_vector(ispec,ispec2D,ipoin, &
- veloc_element, &
+ subroutine wmo_get_vel_vector(ispec,ispec2D, &
+ ipoin,iglob, &
+ val_element, &
narraydim)
! put into this separate routine to make compilation faster
use specfem_par,only: NDIM,ibool
+ use specfem_par_elastic,only: displ,veloc,ispec_is_elastic
+ use specfem_par_acoustic,only: ispec_is_acoustic
use specfem_par_movie
implicit none
- integer :: ispec,ispec2D,ipoin,narraydim
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: &
- veloc_element
+ integer,intent(in) :: ispec,ispec2D,ipoin,iglob
+ integer,intent(in) :: narraydim
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: val_element
! local parameters
- integer :: i,j,k,iglob
+ integer :: i,j,k
logical :: is_done
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc_element(1,i,j,k)
- store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc_element(2,i,j,k)
- store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc_element(3,i,j,k)
- is_done = .true.
- return
- endif
+ ! elastic displacement/velocity
+ if( ispec_is_elastic(ispec) ) then
+ if(SAVE_DISPLACEMENT) then
+ ! velocity x,y,z-components
+ store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin) = displ(1,iglob)
+ store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin) = displ(2,iglob)
+ store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin) = displ(3,iglob)
+ else
+ ! velocity x,y,z-components
+ store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc(1,iglob)
+ store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc(2,iglob)
+ store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc(3,iglob)
+ endif
+ endif
+
+ ! acoustic pressure potential
+ if( ispec_is_acoustic(ispec) ) then
+ is_done = .false.
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ if( iglob == ibool(i,j,k,ispec) ) then
+ store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin) = val_element(1,i,j,k)
+ store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin) = val_element(2,i,j,k)
+ store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin) = val_element(3,i,j,k)
+ is_done = .true.
+ return
+ endif
+ enddo
enddo
enddo
- enddo
+ endif
end subroutine wmo_get_vel_vector
@@ -621,25 +630,12 @@
j = free_surface_ijk(2,igll,iface)
k = free_surface_ijk(3,igll,iface)
iglob = ibool(i,j,k,ispec)
- ! elastic displacement/velocity
- if( ispec_is_elastic(ispec) ) then
- if(SAVE_DISPLACEMENT) then
- store_val_ux_external_mesh(ipoin) = displ(1,iglob)
- store_val_uy_external_mesh(ipoin) = displ(2,iglob)
- store_val_uz_external_mesh(ipoin) = displ(3,iglob)
- else
- store_val_ux_external_mesh(ipoin) = veloc(1,iglob)
- store_val_uy_external_mesh(ipoin) = veloc(2,iglob)
- store_val_uz_external_mesh(ipoin) = veloc(3,iglob)
- endif
- endif
- ! acoustic pressure potential
- if( ispec_is_acoustic(ispec) ) then
- ! stores values from element
- call wmo_get_val_elem(ispec,ipoin,val_element)
- endif
-
+ ! puts displ/velocity values into storage array
+ call wmo_get_vel_vector(ispec,0, &
+ ipoin,iglob, &
+ val_element, &
+ 0)
enddo
else
imin = minval( free_surface_ijk(1,:,iface) )
@@ -659,25 +655,11 @@
iglob = ibool(iorderi(iloc),iorderj(iloc),kmin,ispec)
endif
- ! elastic displacement/velocity
- if( ispec_is_elastic(ispec) ) then
- if(SAVE_DISPLACEMENT) then
- store_val_ux_external_mesh(ipoin) = displ(1,iglob)
- store_val_uy_external_mesh(ipoin) = displ(2,iglob)
- store_val_uz_external_mesh(ipoin) = displ(3,iglob)
- else
- store_val_ux_external_mesh(ipoin) = veloc(1,iglob)
- store_val_uy_external_mesh(ipoin) = veloc(2,iglob)
- store_val_uz_external_mesh(ipoin) = veloc(3,iglob)
- endif
- endif
-
- ! acoustic pressure potential
- if( ispec_is_acoustic(ispec) ) then
- ! stores values from element
- call wmo_get_val_elem(ispec,ipoin,val_element)
- endif
-
+ ! puts displ/velocity values into storage array
+ call wmo_get_vel_vector(ispec,0, &
+ ipoin,iglob, &
+ val_element, &
+ 0)
enddo ! iloc
endif
enddo ! iface
@@ -751,42 +733,6 @@
!================================================================
- subroutine wmo_get_val_elem(ispec,ipoin,val_element)
-
- ! put into this separate routine to make compilation faster
-
- use specfem_par,only: NDIM,ibool
- use specfem_par_movie
- implicit none
-
- integer :: ispec,ipoin
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: &
- val_element
-
- ! local parameters
- integer :: i,j,k,iglob
- logical :: is_done
-
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- store_val_ux_external_mesh(ipoin) = val_element(1,i,j,k)
- store_val_uy_external_mesh(ipoin) = val_element(2,i,j,k)
- store_val_uz_external_mesh(ipoin) = val_element(3,i,j,k)
- is_done = .true.
- return
- endif
- enddo
- enddo
- enddo
-
- end subroutine wmo_get_val_elem
-
-!=====================================================================
-
subroutine wmo_create_shakemap_o()
! outputs shakemap file
@@ -1035,6 +981,9 @@
character(len=3) :: channel
character(len=1) :: compx,compy,compz
+ !real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: tmpdata
+ !integer :: i,j,k,iglob
+
! gets component characters: X/Y/Z or E/N/Z
call write_channel_name(1,channel)
compx(1:1) = channel(3:3) ! either X or E
@@ -1185,6 +1134,57 @@
!write(27) velocity_movie
!close(27)
+ ! norms
+ !allocate( tmpdata(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ !if( ier /= 0 ) stop 'error allocating tmpdata arrays for movie output'
+ !
+ !if( ELASTIC_SIMULATION ) then
+ ! ! norm of displacement
+ ! do ispec=1,NSPEC_AB
+ ! do k=1,NGLLZ
+ ! do j=1,NGLLY
+ ! do i=1,NGLLX
+ ! iglob = ibool(i,j,k,ispec)
+ ! tmpdata(i,j,k,ispec) = sqrt( displ(1,iglob)**2 + displ(2,iglob)**2 + displ(3,iglob)**2 )
+ ! enddo
+ ! enddo
+ ! enddo
+ ! enddo
+ ! write(outputname,"('/proc',i6.6,'_displ_norm_it',i6.6,'.bin')") myrank,it
+ ! open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ ! if( ier /= 0 ) stop 'error opening file movie output velocity z'
+ ! write(27) tmpdata
+ ! close(27)
+ !
+ ! ! norm of acceleration
+ ! do ispec=1,NSPEC_AB
+ ! do k=1,NGLLZ
+ ! do j=1,NGLLY
+ ! do i=1,NGLLX
+ ! iglob = ibool(i,j,k,ispec)
+ ! tmpdata(i,j,k,ispec) = sqrt( accel(1,iglob)**2 + accel(2,iglob)**2 + accel(3,iglob)**2 )
+ ! enddo
+ ! enddo
+ ! enddo
+ ! enddo
+ ! write(outputname,"('/proc',i6.6,'_accel_norm_it',i6.6,'.bin')") myrank,it
+ ! open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ ! if( ier /= 0 ) stop 'error opening file movie output velocity z'
+ ! write(27) tmpdata
+ ! close(27)
+ !endif
+ !
+ !! norm of velocity
+ !tmpdata = sqrt( velocity_x**2 + velocity_y**2 + velocity_z**2)
+ !
+ !write(outputname,"('/proc',i6.6,'_velocity_norm_it',i6.6,'.bin')") myrank,it
+ !open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ !if( ier /= 0 ) stop 'error opening file movie output velocity z'
+ !write(27) tmpdata
+ !close(27)
+ !
+ !deallocate(tmpdata)
+
endif
end subroutine wmo_movie_volume_output
Added: seismo/3D/SPECFEM3D/trunk/utils/locate_partition.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/locate_partition.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/locate_partition.f90 2011-11-29 15:09:45 UTC (rev 19248)
@@ -0,0 +1,220 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 2 . 0
+! ---------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA and University of Pau / CNRS / INRIA
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! April 2011
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! utility to locate partition which is closest to given point location
+!
+! compile with:
+! > gfortran -I src/shared -o bin/xlocate_partition utils/locate_partition.f90
+! or
+! > ifort -assume byterecl -I src/shared -o bin/xlocate_partition utils/locate_partition.f90
+!
+! run with:
+! > ./bin/xlocate_partition 70000.0 11000.0 -3000.0 ./in_out_files/DATABASES_MPI/
+
+ program locate_partition
+
+! works for external, unregular meshes
+
+ implicit none
+
+ include 'constants.h'
+
+ ! mesh coordinates
+ real(kind=CUSTOM_REAL),dimension(:),allocatable :: xstore, ystore, zstore
+ integer, dimension(:,:,:,:),allocatable :: ibool
+ integer :: NSPEC_AB, NGLOB_AB
+
+ integer :: i,ios,ier
+ integer :: iproc
+ character(len=256) :: arg(4)
+ character(len=256) :: LOCAL_PATH
+ character(len=256) :: prname_lp
+
+ real(kind=CUSTOM_REAL) :: x_found,y_found,z_found,distance
+ real(kind=CUSTOM_REAL) :: target_x,target_y,target_z
+ real(kind=CUSTOM_REAL) :: total_x,total_y,total_z
+ real(kind=CUSTOM_REAL) :: total_distance
+ integer :: total_partition
+
+! checks given arguments
+ print *
+ print *,'locate partition'
+ print *,'----------------------------'
+
+ do i = 1, 4
+ call getarg(i,arg(i))
+ if (i <= 4 .and. trim(arg(i)) == '') then
+ print *, 'Usage: '
+ print *, ' xlocate_partition x y z Databases_directory'
+ stop ' Reenter command line options'
+ endif
+ enddo
+ read(arg(1),*) target_x
+ read(arg(2),*) target_y
+ read(arg(3),*) target_z
+ LOCAL_PATH = arg(4)
+
+ print *,'search location: '
+ print *,' x = ',target_x
+ print *,' y = ',target_y
+ print *,' z = ',target_z
+ print *,'in directory: ',trim(LOCAL_PATH)
+ print *,'----------------------------'
+ print *
+
+ ! loops over slices (process partitions)
+ total_distance = HUGEVAL
+ total_partition = -1
+ total_x = 0.0
+ total_y = 0.0
+ total_z = 0.0
+ iproc = -1
+ do while( iproc < 10000000 )
+ ! starts with 0
+ iproc = iproc + 1
+
+ ! gets number of elements and global points for this partition
+ write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'/proc',iproc,'_'
+ open(unit=27,file=prname_lp(1:len_trim(prname_lp))//'external_mesh.bin',&
+ status='old',action='read',form='unformatted',iostat=ios)
+ if( ios /= 0 ) exit
+
+ read(27,iostat=ier) NSPEC_AB
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+ read(27,iostat=ier) NGLOB_AB
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+
+ ! ibool file
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
+ read(27,iostat=ier) ibool
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+
+ ! global point arrays
+ allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xstore etc.'
+ read(27,iostat=ier) xstore
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+ read(27,iostat=ier) ystore
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+ read(27,iostat=ier) zstore
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+ close(27)
+
+ print*,'partition: ',iproc
+ print*,' min/max x = ',minval(xstore),maxval(xstore)
+ print*,' min/max y = ',minval(ystore),maxval(ystore)
+ print*,' min/max z = ',minval(zstore),maxval(zstore)
+ print*
+
+ ! gets distance to target location
+ call get_closest_point(target_x,target_y,target_z, &
+ NGLOB_AB,NSPEC_AB,xstore,ystore,zstore,ibool, &
+ distance,x_found,y_found,z_found)
+
+ if( distance < total_distance ) then
+ total_distance = distance
+ total_partition = iproc
+ total_x = x_found
+ total_y = y_found
+ total_z = z_found
+ endif
+
+ ! cleans up memory allocations
+ deallocate(ibool,xstore,ystore,zstore)
+
+ enddo ! all slices for points
+
+ ! checks
+ if (total_partition < 0 ) then
+ print*,'Error: partition not found among ',iproc,'partitions searched'
+ stop 'Error: partition not found'
+ endif
+
+ ! output
+ print*,'number of partitions searched: ',iproc
+ print*
+ print*,'closest grid point location found:'
+ print*,' x = ',total_x
+ print*,' y = ',total_y
+ print*,' z = ',total_z
+ print*,' distance to search location = ',sqrt(total_distance)
+ print*,'closest partition: '
+ print*,' partition = ',total_partition
+ print*
+
+ end program locate_partition
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine get_closest_point(target_x,target_y,target_z, &
+ NGLOB_AB,NSPEC_AB,xstore,ystore,zstore,ibool, &
+ distance,x_found,y_found,z_found)
+
+ implicit none
+ include 'constants.h'
+
+ real(kind=CUSTOM_REAL),intent(in) :: target_x,target_y,target_z
+
+ integer,intent(in) :: NSPEC_AB,NGLOB_AB
+ real(kind=CUSTOM_REAL),dimension(NGLOB_AB),intent(in) :: xstore, ystore, zstore
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: ibool
+
+ real(kind=CUSTOM_REAL),intent(out) :: distance
+ real(kind=CUSTOM_REAL),intent(out) :: x_found,y_found,z_found
+
+ ! local parameters
+ integer :: i,j,k,ispec,iglob
+ real(kind=CUSTOM_REAL) :: dist
+
+ distance = HUGEVAL
+ x_found = 0.0
+ y_found = 0.0
+ z_found = 0.0
+
+ do ispec=1,NSPEC_AB
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dist = (target_x - xstore(iglob))*(target_x - xstore(iglob)) &
+ + (target_y - ystore(iglob))*(target_y - ystore(iglob)) &
+ + (target_z - zstore(iglob))*(target_z - zstore(iglob))
+
+ if( dist < distance ) then
+ distance = dist
+ x_found = xstore(iglob)
+ y_found = ystore(iglob)
+ z_found = zstore(iglob)
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+
+ end subroutine
More information about the CIG-COMMITS
mailing list