from rsf.proj import *
from rsf.recipes.m1d import *
nn = 8
mode = 'em'
spmin = 0.0
spmax = 5.0
recmin = 0.6
recmax = 1.0
model = {
'X' : 5,
'T' : 1.0,
'dx': 0.005,
'dt': 0.0005,
'SelT': 0.8,
'snpintvl': 1.0,
'size' : 12,
'frqcut' : 0.6,
'pml' : 240,
'vel' : '2.1+0.1*x1*x1',
'dvel' : '0.2*x1',
'den' : '2'
}
srp = {
'bgn' : 0.15,
'frq' : 12.5,
'srcmms' : 'y',
'inject': 'y',
'slx' : 0.5,
'gdep' : 2.5
}
def modeling(errpar, ipar, mode, idn):
dx = ipar['dx']
dt = ipar['dt']
if mode == 'et':
_pf = 't'
var = ipar['dt']
elif mode == 'em':
_pf = 'm'
var = ipar['dx']
elif mode == 'ea':
_pf = 'a'
var = ipar['dt']+ipar['dx']
vel = '%svel%d' %(_pf, idn)
dvel = '%sdvel%d' %(_pf, idn)
den = '%sden%d' %(_pf, idn)
velhf = vel+'hf'
dvelhf= dvel+'hf'
denhf = den+'hf'
vel_pml = vel+'_pml'
den_pml = den+'_pml'
buildmodel(ipar, den, vel, dvel, ipar['den'], ipar['vel'], ipar['dvel'])
srcp = '%ssrcp%d' %(_pf, idn)
srcd = '%ssrcd%d' %(_pf, idn)
presrc = '%spresrc%d' %(_pf, idn)
velsrc = '%svelsrc%d' %(_pf, idn)
preinit = '%spreinit%d' %(_pf, idn)
velinit = '%svelinit%d' %(_pf, idn)
mms = '%smms%d' %(_pf, idn)
if ipar['inject'] == 'n' :
aslt = '%saslt%d' %(_pf, idn)
analyticslt(aslt, ipar, float(ipar['vel']), _pf, idn)
mmsfiles = {}
if ipar['srcmms'] == 'y' :
buildsrcd(ipar, srcd, _pf, idn)
buildmms(ipar, mms, presrc, velsrc, preinit, velinit, vel, dvel, den, velhf, dvelhf, denhf)
mmsfiles = {'presrc' : presrc, 'velsrc' :velsrc,
'preinit': preinit, 'velinit':velinit}
mgraph(mms,'MMSSLT-%s-%g'%(mode,var))
mmssnap = '%smmssnap%d' %(_pf, idn)
Flow(mmssnap, mms, 'window n2=1 f2=%(snt)d ' %ipar)
graph(spmin, spmax, mmssnap,"MMS-Snap-%s-%g" %(mode,var))
mmsrec = '%smmsrec%d' %(_pf, idn)
Flow(mmsrec, mms, 'window n1=1 f1=%(gp)d' %ipar)
graph(recmin, recmax, mmsrec,"MMS-Rec-%s-%g" %(mode,var))
buildsrcp(ipar, srcp)
ic = '%sic%d' %(_pf, idn)
buildic(ipar, ic)
lrw = '%slrwav%d' %(_pf, idn)
lrr = '%slrrec%d' %(_pf, idn)
lrs = '%slrsnap%d' %(_pf, idn)
lrmodel(lrw, lrr, srcp, ic, vel, den, mmsfiles, ipar, _pf, idn)
mgraph(lrw,'LR-%s-%g'%(mode,var))
Flow(lrs, lrw, 'window n2=1 f2=%(snt)d ' %ipar)
graph(13.5, 14.5, lrs,"LR-Rec-%s-%g" %(mode,var))
graph(recmin, recmax, lrr,"LR-Rec-%s-%g" %(mode,var))
lfdw = '%slfdwav%d' %(_pf, idn)
lfdr = '%slfdrec%d' %(_pf, idn)
lfds = '%slfdsnap%d' %(_pf, idn)
lfdmodel(lfdw, lfdr, srcp, ic , vel_pml, den_pml, mmsfiles, ipar, _pf, idn)
mgraph(lfdw,'LFD-%s-%g'%(mode,var))
Flow(lfds, lfdw, 'window n2=1 f2=%(snt)d ' %ipar)
graph(spmin, spmax, lfds,"LFD-Rec-%s-%g" %(mode,var))
graph(recmin, recmax, lfdr,"LFD-Rec-%s-%g" %(mode,var))
fdw = '%sfdwav%d' %(_pf, idn)
fdr = '%sfdrec%d' %(_pf, idn)
fds = '%sfdsnap%d' %(_pf, idn)
fdmodel(fdw, fdr, srcp, ic, vel_pml, den_pml, mmsfiles, ipar)
mgraph(fdw,'FD-%s-%g'%(mode,var))
Flow(fds, fdw, 'window n2=1 f2=%(snt)d ' %ipar)
graph(spmin, spmax, fds,"FD-Rec-%s-%g" %(mode,var))
graph(recmin, recmax, fdr,"FD-Rec-%s-%g" %(mode,var))
snap = '%ssnap%d' %(_pf, idn)
if ipar['srcmms'] == 'y' and ipar['inject']=='y' :
snaplist=[mmssnap, lrs, lfds, fds]
elif ipar['inject']=='n':
snaplist=[aslt, lrs, lfds, fds]
else :
snaplist=[lrs, lfds, fds]
snapplot(snap, snaplist, spmin,spmax )
recplot = '%srec%d' %(_pf, idn)
if ipar['srcmms'] == 'y' and ipar['inject']=='y' :
reclist=[mmsrec, lrr, lfdr, fdr]
elif ipar['inject']=='n':
reclist=[lrr, lfdr, fdr]
else :
reclist=[lrr, lfdr, fdr]
cat(recplot+'_',reclist,2)
Flow(recplot,recplot+'_','put label2="Amplitude" unit2=""')
graph(recmin,recmax,recplot,recplot)
if ipar['srcmms'] == 'n' and ipar['inject']=='y':
lrwav0 = '%slrsnap0' %_pf
bmk = '%sbanckmark%d' %(_pf, idn)
if mode == 'et':
sample(bmk,lrwav0, idn, 2)
elif mode == 'em':
sample(bmk,lrwav0, idn, 1)
elif mode == 'ea':
sample(bmk,lrwav0, idn, 3)
else:
print "Error of mode!"
emit
elif ipar['inject']=='n' :
bmk = aslt
elif ipar['srcmms'] == 'y' and ipar['inject']=='y' :
bmk = mmssnap
lre = '%slrerr%d' %(_pf, idn)
err(lre, lrs, bmk)
lfde = '%slfderr%d' %(_pf, idn)
err(lfde, lfds, bmk)
fde = '%sfderr%d' %(_pf, idn)
err(fde, fds, bmk)
if errpar == {}:
errpar['%slr' %_pf] = [lre]
errpar['%slfd' %_pf] = [lfde]
errpar['%sfd' %_pf] = [fde]
else:
errpar['%slr' %_pf].append(lre)
errpar['%slfd' %_pf].append(lfde)
errpar['%sfd' %_pf].append(fde)
def sample(fout, fin, smp, axis):
jj = 2**smp
if axis == '1' or axis == 1 :
Flow(fout, fin,
'''
window j1=%d
''' %jj)
elif axis == '2' or axis == 2:
Flow(fout, fin,
'''
window j2=%d
''' %jj)
elif axis == '3' or axis == 3:
Flow(fout, fin,
'''
window j1=%d j2=%d
''' %(jj,jj))
def err(fout, fin, banchmark):
Flow(fout, [fin, banchmark],
'''
add ${SOURCES[1]} scale=1,-1 |
stack axis=1 rms=y |
stack axis=1 rms=y |
math output="(input)"
''')
def cat(fout, flist, axis):
ni = len(flist)
src = ""
for ii in range(1,ni):
src += "${SOURCES[%d]} " %ii
Flow(fout, flist, 'cat %s axis=%d ' %(src,axis))
def errcur(errc, err, par, mode, nn):
dx = par['dx']
dt = par['dt']
frq = par['frq']
if mode == 'et':
_pf = 't'
axis = '%sax' %_pf
Flow(axis, None,
'''
spike n1=%d o1=0.0 d1=1.0 |
put label1="dt" unit1="ms" |
math output="(%g*(x1+1))*1000"
'''%(nn, dt))
elif mode == 'em':
_pf = 'm'
axis = '%sax' %_pf
Flow(axis, None,
'''
spike n1=%d o1=0.0 d1=1.0 |
put label1="dx" unit1="m" |
math output="(%g*(x1+1))*1000"
'''%(nn, dx))
elif mode == 'ea':
_pf = 'a'
axis = '%sax' %_pf
Flow(axis, None,
'''
spike n1=%d o1=0.0 d1=1.0 |
put label1="dx" unit1="m" |
math output="(%g*(2.0^x1))*1000"
'''%(nn, dx))
else :
axis = ''
Flow(errc,[axis,err],'cmplx ${SOURCES[1]} | put label2="Error(RMS)"')
def mgraph(fin, title):
Plot(fin,
'''
put label1 = "depth" unit1="m" label2="Amplitude" unit2="" |
window j2=50 squeeze=n |
transp plane=23|
graph grid1=y grid2=y max2=1.0 min2=-1.0 title=%s
'''%(str(title)))
def graph(min1, max1, fin, title):
Plot(fin,
'''
graph grid1=n grid2=n g2num=0.25 g1num=0.1
wanttitle=n
plotcol=7,6,3,5 dash=0,1,2,3,5
min1=%g max1=%g min2=-0.6 max2=1.2
screenratio=0.618 screenht=8 plotfat=4
'''%(float(min1),float(max1)))
def snapplot(fout, flist, min2, max2):
cat(fout,flist,2)
graph(min2,max2,fout,fout)
par0 = setpar(model, srp)
par = par0
errlist = {}
dx0 = model['dx']
dt0 = model['dt']
for ii in range(0, nn):
modeling(errlist, par, mode, ii)
if mode == 'et' :
model['dt'] = model['dt']+dt0
elif mode == 'em':
model['dx'] = model['dx']+dx0
elif mode == 'ea':
model['dt'] = model['dt']*2
model['dx'] = model['dx']*2
par = setpar(model, srp)
lcur = []
keys = errlist.keys()
for key in sorted(keys,reverse=True):
cat(key, errlist[key][0:nn], 1)
Flow(key+'s',key,'sfsmooth rect1=2')
errcur(key+'c', key+'s', par0, mode, nn)
lcur += [key+'c']
cat(mode,lcur, 2)
Result(mode,
'''
graph min2=-0.001 max2=0.03 max1=40.5
wanttitle=n grid1=n grid2=n g2num=0.01
screenratio=0.6 screenht=8
labelsz=8
plotfat=4 plotcol=6,3,5 dash=3,2,0
''')
print errlist.keys()
End() |
![[pdf]](paper_html/icons/pdf.png){kind=icon}