forked from SPECFEM/specfem3d
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlocate_partition.f90
230 lines (192 loc) · 7.76 KB
/
locate_partition.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
!=====================================================================
!
! S p e c f e m 3 D V e r s i o n 3 . 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 3 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 one of these (use your default):
! > gfortran -I src/shared -o bin/xlocate_partition utils/locate_partition.f90
! > ifort -assume byterecl -I src/shared -o bin/xlocate_partition utils/locate_partition.f90
! > pgf90 -I src/shared -o bin/xlocate_partition utils/locate_partition.f90
!
! specify a target (x,y) may be in UTM, not lon-lat, then run with:
! > ./bin/xlocate_partition 70000.0 11000.0 -3000.0 ./OUTPUT_FILES/DATABASES_MPI/
!
! this will generate the output file OUTPUT_FILES/DATABASES_MPI/partition_bounds.dat
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 get_command_argument(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 *
! open a text file to list the maximal bounds of each partition
open(11,file=trim(LOCAL_PATH)//'partition_bounds.dat',status='unknown')
! 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'
read(27,iostat=ier) ios ! skip dummy
! 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 *
write(11,'(i10,6e18.6)') iproc,minval(xstore),maxval(xstore),minval(ystore),maxval(ystore),minval(zstore),maxval(zstore)
! 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
close(11)
! 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