Skip to content

Commit 8cc8580

Browse files
author
Diego Melgar
committed
Added MPI support for synthetics
1 parent 7e5eaa5 commit 8cc8580

2 files changed

Lines changed: 407 additions & 79 deletions

File tree

‎src/python/mudpy/parallel.py‎

Lines changed: 244 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
'''
44

55

6-
def run_paralell_green(home,project_name,station_file,model_name,dt,NFFT,static,dk,pmin,pmax,kmax,tsunami,rank,size):
6+
def run_parallel_green(home,project_name,station_file,model_name,dt,NFFT,static,dk,pmin,pmax,kmax,tsunami,rank,size):
77
'''
88
Compute GFs using Zhu & Rivera code for a given velocity model, source depth
99
and station file. This function will make an external system call to fk.pl
@@ -87,18 +87,247 @@ def run_paralell_green(home,project_name,station_file,model_name,dt,NFFT,static,
8787
p=subprocess.Popen(command,stdout=open(write_file,'w'),stderr=subprocess.PIPE)
8888
p.communicate()
8989

90+
91+
def run_parallel_synthetics(home,project_name,station_file,model_name,integrate,static,tsunami,time_epi,beta,rank,size):
92+
'''
93+
Use green functions and compute synthetics at stations for a single source
94+
and multiple stations. This code makes an external system call to syn.c first it
95+
will make the external call for the strike-slip component then a second externall
96+
call will be made for the dip-slip component. The unit amount of moment is 1e15
97+
which corresponds to Mw=3.9333...
98+
99+
IN:
100+
source: 1-row numpy array containig informaiton aboutt he source, lat, lon, depth, etc...
101+
station_file: File name with the station coordinates
102+
green_path: Directopry where GFs are stored
103+
model_file: File containing the Earth velocity structure
104+
integrate: =0 if youw ant velocity waveforms, =1 if you want displacements
105+
static: =0 if computing full waveforms, =1 if computing only the static field
106+
subfault: String indicating the subfault being worked on
107+
coord_type: =0 if problem is in cartesian coordinates, =1 if problem is in lat/lon
108+
109+
OUT:
110+
log: Sysytem standard output and standard error for log
111+
'''
112+
113+
import os
114+
import subprocess
115+
from mudpy.forward import get_mu
116+
from string import rjust
117+
from numpy import array,genfromtxt,loadtxt,savetxt,log10
118+
from obspy import read
119+
from shlex import split
120+
from mudpy.green import src2sta,rt2ne,origin_time
121+
from progressbar import ProgressBar
122+
123+
#What parameters are we using?
124+
if rank==0:
125+
out='''Running all processes with:
126+
home = %s
127+
project_name = %s
128+
station_file = %s
129+
model_name = %s
130+
integrate = %s
131+
static = %s
132+
tsunami = %s
133+
time_epi = %s
134+
beta = %d
135+
''' %(home,project_name,station_file,model_name,str(integrate),str(static),str(tsunami),str(time_epi),beta)
136+
print out
137+
#Read your corresponding source file
138+
mpi_source=genfromtxt(home+project_name+'/data/model_info/mpi_source.'+str(rank)+'.fault')
139+
#Constant parameters
140+
rakeDS=90+beta #90 is thrust, -90 is normal
141+
rakeSS=0+beta #0 is left lateral, 180 is right lateral
142+
tb=50 #Number of samples before first arrival
143+
#Load structure
144+
model_file=home+project_name+'/structure/'+model_name
145+
structure=loadtxt(model_file,ndmin=2)
146+
for ksource in range(len(mpi_source)):
147+
source=mpi_source[ksource,:]
148+
#Parse the soruce information
149+
num=rjust(str(int(source[0])),4,'0')
150+
zs=source[3]
151+
strike=source[4]
152+
dip=source[5]
153+
rise=source[6]
154+
duration=source[7]
155+
ss_length=source[8]
156+
ds_length=source[9]
157+
strdepth='%.4f' % zs
158+
subfault=rjust(str(int(source[0])),4,'0')
159+
if static==0 and tsunami==0: #Where to save dynamic waveforms
160+
green_path=home+project_name+'/GFs/dynamic/'+model_name+"_"+strdepth+".sub"+subfault+"/"
161+
if static==0 and tsunami==1: #Where to save dynamic waveforms
162+
green_path=home+project_name+'/GFs/tsunami/'+model_name+"_"+strdepth+".sub"+subfault+"/"
163+
if static==1: #Where to save dynamic waveforms
164+
green_path=home+project_name+'/GFs/static/'
165+
staname=genfromtxt(station_file,dtype="S6",usecols=0)
166+
if staname.shape==(): #Single staiton file
167+
staname=array([staname])
168+
#Compute distances and azmuths
169+
d,az=src2sta(station_file,source)
170+
#Get moment corresponding to 1 meter of slip on subfault
171+
mu=get_mu(structure,zs)
172+
Mo=mu*ss_length*ds_length*1
173+
Mw=(2./3)*(log10(Mo)-9.1)
174+
#Move to output folder
175+
os.chdir(green_path)
176+
#Instantiate progress bar
177+
pbar = ProgressBar()
178+
print 'Processor '+str(rank)+' is working on subfault '+str(int(source[0]))+' and '+str(len(d))+' stations '
179+
for k in range(len(d)):
180+
if static==0: #Compute full waveforms
181+
diststr='%.3f' % d[k] #Need current distance in string form for external call
182+
#Form the strings to be used for the system calls according to suer desired options
183+
if integrate==1: #Make displ.
184+
#First Stike-Slip GFs
185+
commandSS="syn -I -M"+str(Mw)+"/"+str(strike)+"/"+str(dip)+"/"+str(rakeSS)+" -D"+str(duration)+ \
186+
"/"+str(rise)+" -A"+str(az[k])+" -O"+staname[k]+".subfault"+num+".SS.disp.x -G"+green_path+diststr+".grn.0"
187+
commandSS=split(commandSS) #Split string into lexical components for system call
188+
#Now dip slip
189+
commandDS="syn -I -M"+str(Mw)+"/"+str(strike)+"/"+str(dip)+"/"+str(rakeDS)+" -D"+str(duration)+ \
190+
"/"+str(rise)+" -A"+str(az[k])+" -O"+staname[k]+".subfault"+num+".DS.disp.x -G"+green_path+diststr+".grn.0"
191+
commandDS=split(commandDS)
192+
else: #Make vel.
193+
#First Stike-Slip GFs
194+
commandSS="syn -M"+str(Mw)+"/"+str(strike)+"/"+str(dip)+"/"+str(rakeSS)+" -D"+str(duration)+ \
195+
"/"+str(rise)+" -A"+str(az[k])+" -O"+staname[k]+".subfault"+num+".SS.vel.x -G"+green_path+diststr+".grn.0"
196+
commandSS=split(commandSS)
197+
#Now dip slip
198+
commandDS="syn -M"+str(Mw)+"/"+str(strike)+"/"+str(dip)+"/"+str(rakeDS)+" -D"+str(duration)+ \
199+
"/"+str(rise)+" -A"+str(az[k])+" -O"+staname[k]+".subfault"+num+".DS.vel.x -G"+green_path+diststr+".grn.0"
200+
commandDS=split(commandDS)
201+
#Run the strike- and dip-slip commands (make system calls)
202+
p=subprocess.Popen(commandSS)
203+
p.communicate()
204+
p=subprocess.Popen(commandDS)
205+
p.communicate()
206+
#Result is in RTZ system (+Z is down) rotate to NEZ with +Z up and scale to m or m/s
207+
if integrate==1: #'tis displacememnt
208+
#Strike slip
209+
r=read(staname[k]+".subfault"+num+'.SS.disp.r')
210+
t=read(staname[k]+".subfault"+num+'.SS.disp.t')
211+
z=read(staname[k]+".subfault"+num+'.SS.disp.z')
212+
ntemp,etemp=rt2ne(r[0].data,t[0].data,az[k])
213+
#Scale to m and overwrite with rotated waveforms
214+
n=r.copy()
215+
n[0].data=ntemp/100
216+
e=t.copy()
217+
e[0].data=etemp/100
218+
z[0].data=z[0].data/100
219+
n=origin_time(n,time_epi,tb)
220+
e=origin_time(e,time_epi,tb)
221+
z=origin_time(z,time_epi,tb)
222+
n.write(staname[k]+".subfault"+num+'.SS.disp.n',format='SAC')
223+
e.write(staname[k]+".subfault"+num+'.SS.disp.e',format='SAC')
224+
z.write(staname[k]+".subfault"+num+'.SS.disp.z',format='SAC')
225+
#Dip Slip
226+
r=read(staname[k]+".subfault"+num+'.DS.disp.r')
227+
t=read(staname[k]+".subfault"+num+'.DS.disp.t')
228+
z=read(staname[k]+".subfault"+num+'.DS.disp.z')
229+
ntemp,etemp=rt2ne(r[0].data,t[0].data,az[k])
230+
n=r.copy()
231+
n[0].data=ntemp/100
232+
e=t.copy()
233+
e[0].data=etemp/100
234+
z[0].data=z[0].data/100
235+
n=origin_time(n,time_epi,tb)
236+
e=origin_time(e,time_epi,tb)
237+
z=origin_time(z,time_epi,tb)
238+
n.write(staname[k]+".subfault"+num+'.DS.disp.n',format='SAC')
239+
e.write(staname[k]+".subfault"+num+'.DS.disp.e',format='SAC')
240+
z.write(staname[k]+".subfault"+num+'.DS.disp.z',format='SAC')
241+
else: #Waveforms are velocity, as before, rotate from RT-Z to NE+Z and scale to m/s
242+
#Strike slip
243+
r=read(staname[k]+".subfault"+num+'.SS.vel.r')
244+
t=read(staname[k]+".subfault"+num+'.SS.vel.t')
245+
z=read(staname[k]+".subfault"+num+'.SS.vel.z')
246+
ntemp,etemp=rt2ne(r[0].data,t[0].data,az[k])
247+
n=r.copy()
248+
n[0].data=ntemp/100
249+
e=t.copy()
250+
e[0].data=etemp/100
251+
z[0].data=z[0].data/100
252+
n=origin_time(n,time_epi,tb)
253+
e=origin_time(e,time_epi,tb)
254+
z=origin_time(z,time_epi,tb)
255+
n.write(staname[k]+".subfault"+num+'.SS.vel.n',format='SAC')
256+
e.write(staname[k]+".subfault"+num+'.SS.vel.e',format='SAC')
257+
z.write(staname[k]+".subfault"+num+'.SS.vel.z',format='SAC')
258+
#Dip Slip
259+
r=read(staname[k]+".subfault"+num+'.DS.vel.r')
260+
t=read(staname[k]+".subfault"+num+'.DS.vel.t')
261+
z=read(staname[k]+".subfault"+num+'.DS.vel.z')
262+
ntemp,etemp=rt2ne(r[0].data,t[0].data,az[k])
263+
n=r.copy()
264+
n[0].data=ntemp/100
265+
e=t.copy()
266+
e[0].data=etemp/100
267+
z[0].data=z[0].data/100
268+
n=origin_time(n,time_epi,tb)
269+
e=origin_time(e,time_epi,tb)
270+
z=origin_time(z,time_epi,tb)
271+
n.write(staname[k]+".subfault"+num+'.DS.vel.n',format='SAC')
272+
e.write(staname[k]+".subfault"+num+'.DS.vel.e',format='SAC')
273+
z.write(staname[k]+".subfault"+num+'.DS.vel.z',format='SAC')
274+
else: #Compute static synthetics
275+
diststr='%.1f' % d[k] #Need current distance in string form for external call
276+
green_file=model_name+".static."+strdepth+".sub"+subfault #Output dir
277+
statics=loadtxt(green_file) #Load GFs
278+
#Print static GFs into a pipe and pass into synthetics command
279+
try:
280+
temp_pipe=statics[k,:]
281+
except:
282+
temp_pipe=statics
283+
inpipe=''
284+
for j in range(len(temp_pipe)):
285+
inpipe=inpipe+' %.6e' % temp_pipe[j]
286+
#Form command for external call
287+
commandDS="syn -M"+str(Mw)+"/"+str(strike)+"/"+str(dip)+"/"+str(rakeDS)+\
288+
" -A"+str(az[k])+" -P"
289+
commandSS="syn -M"+str(Mw)+"/"+str(strike)+"/"+str(dip)+"/"+str(rakeSS)+\
290+
" -A"+str(az[k])+" -P"
291+
commandSS=split(commandSS) #Lexical split
292+
commandDS=split(commandDS)
293+
#Make system calls, one for DS, one for SS, and save log
294+
ps=subprocess.Popen(['printf',inpipe],stdout=subprocess.PIPE) #This is the statics pipe, pint stdout to syn's stdin
295+
p=subprocess.Popen(commandSS,stdin=ps.stdout,stdout=open(staname[k]+'.subfault'+num+'.SS.static.rtz','w'))
296+
p.communicate()
297+
ps=subprocess.Popen(['printf',inpipe],stdout=subprocess.PIPE) #This is the statics pipe, pint stdout to syn's stdin
298+
p=subprocess.Popen(commandDS,stdin=ps.stdout,stdout=open(staname[k]+'.subfault'+num+'.DS.static.rtz','w'))
299+
p.communicate()
300+
#Rotate radial/transverse to East/North, correct vertical and scale to m
301+
statics=loadtxt(staname[k]+'.subfault'+num+'.SS.static.rtz')
302+
u=statics[2]/100
303+
r=statics[3]/100
304+
t=statics[4]/100
305+
ntemp,etemp=rt2ne(array([r,r]),array([t,t]),az[k])
306+
n=ntemp[0]
307+
e=etemp[0]
308+
savetxt(staname[k]+'.subfault'+num+'.SS.static.neu',(n,e,u,beta))
309+
statics=loadtxt(staname[k]+'.subfault'+num+'.DS.static.rtz')
310+
u=statics[2]/100
311+
r=statics[3]/100
312+
t=statics[4]/100
313+
ntemp,etemp=rt2ne(array([r,r]),array([t,t]),az[k])
314+
n=ntemp[0]
315+
e=etemp[0]
316+
savetxt(staname[k]+'.subfault'+num+'.DS.static.neu',(n,e,u,beta))
317+
90318

91319

92320
#If main entry point
93321
if __name__ == '__main__':
94322
import sys
95323
from mpi4py import MPI
324+
from obspy import UTCDateTime
96325

97326
comm = MPI.COMM_WORLD
98327
rank = comm.Get_rank()
99328
size = comm.Get_size()
100329
# Map command line arguments to function arguments.
101-
if sys.argv[1]=='run_paralell_green':
330+
if sys.argv[1]=='run_parallel_green':
102331
#Parse command line arguments
103332
home=sys.argv[2]
104333
project_name=sys.argv[3]
@@ -112,6 +341,18 @@ def run_paralell_green(home,project_name,station_file,model_name,dt,NFFT,static,
112341
pmax=int(sys.argv[11])
113342
kmax=int(sys.argv[12])
114343
tsunami=sys.argv[13]=='True'
115-
run_paralell_green(home,project_name,station_file,model_name,dt,NFFT,static,dk,pmin,pmax,kmax,tsunami,rank,size)
344+
run_parallel_green(home,project_name,station_file,model_name,dt,NFFT,static,dk,pmin,pmax,kmax,tsunami,rank,size)
345+
elif sys.argv[1]=='run_parallel_synthetics':
346+
#Parse command line arguments
347+
home=sys.argv[2]
348+
project_name=sys.argv[3]
349+
station_file=sys.argv[4]
350+
model_name=sys.argv[5]
351+
integrate=int(sys.argv[6])
352+
static=int(sys.argv[7])
353+
tsunami=sys.argv[8]=='True'
354+
time_epi=UTCDateTime(sys.argv[9])
355+
beta=float(sys.argv[10])
356+
run_parallel_synthetics(home,project_name,station_file,model_name,integrate,static,tsunami,time_epi,beta,rank,size)
116357
else:
117358
print 'ERROR: You''re not allowed to run '+sys.argv[1]+' from the shell or it does not exist'

0 commit comments

Comments
 (0)