# Sconstruct for 3D esg: (1) record trace data; # (2) generate movie; # (3) display discrete total energy of the solution as a function of time # (serial version) # Author: MZ from rsf.proj import* CWPROOT = os.getenv('CWPROOT') sunull = os.path.join(CWPROOT,'bin/sunull') sushw = os.path.join(CWPROOT,'bin/sushw') suwaveform = os.path.join(CWPROOT,'bin/suwaveform') RSFSRC = os.getenv('RSFSRC') esg = os.path.join(RSFSRC,'trip2.1/iwave/esg/main/esg.x') stdmdl = os.path.join(RSFSRC,'trip2.1/iwave/grid/main/standardmodel.x') Flow('cout0.txt',None,'touch $TARGET') ### === SET PARAMETERS import math ## for source fpeak=20 T=0.2 ## for stiffness fields rho=1 cs=1 cp=1.5 mu0=rho*cs**2; k0=rho*cp**2-4*mu0/3; ## for grid nz=75 nx=nz+4 ny=nz-3 dz=cs*1000/fpeak/10 dx=0.98*dz dy=0.97*dz dt=min(dz,dx,dy)/cp/math.sqrt(3)*0.5 ### === CREATE FIELDS Flow('byh', None, ''' makevel n1=%d n2=%d n3=%d d1=%f d2=%f d3=%f v000=1 | sfput dim=3 gdim=3 id1=0 id2=1 id3=2 '''%(nz,nx,ny,dz,dx,dy), stdin=0) Flow('c11h', None, ''' makevel n1=%d n2=%d n3=%d d1=%f d2=%f d3=%f v000=%f | sfput dim=3 gdim=3 id1=0 id2=1 id3=2 '''%(nz,nx,ny,dz,dx,dy,k0+4*mu0/3), stdin=0) Flow('c13h', None, ''' makevel n1=%d n2=%d n3=%d d1=%f d2=%f d3=%f v000=%f | sfput dim=3 gdim=3 id1=0 id2=1 id3=2 '''%(nz,nx,ny,dz,dx,dy,k0-2*mu0/3), stdin=0) Flow('c55h', None, ''' makevel n1=%d n2=%d n3=%d d1=%f d2=%f d3=%f v000=%f | sfput dim=3 gdim=3 id1=0 id2=1 id3=2 '''%(nz,nx,ny,dz,dx,dy,mu0), stdin=0) ### === CREATE RECEIVER TRACE FILE Flow('data.su', None, sunull + ' nt=301 ntr=%d dt=0.002 | '%(nx/2) + sushw + ' key=gx a=%g b=%g | '%(nx/4*dx,dx) + sushw + ' key=gy a=%g | '%(math.ceil(dy*ny/2)) + sushw + ' key=gelev a=%g > ${TARGETS[0]}'%(-math.ceil(dz*nz/3)), stdin=0) ### === CREATE SOURCE t0=1/fpeak nt=math.ceil(T*1000/dt) # IWAVE presumes that a source is located at its receiver position. Flow('src.su', None, suwaveform + ' type=ricker1 fpeak=%g ns=%d | '%(fpeak, nt) + sushw + ' key=delrt a=%g | '%(t0) + sushw + ' key=gx a=%g | '%(math.ceil(dx*nx/2)) + sushw + ' key=gy a=%g | '%(math.ceil(dy*ny/2)) + sushw + ' key=gelev a=%g > ${TARGETS[0]}'%(-math.ceil(dz*nz/2)), stdin=0) ### === START ## (1) create data_vx profile Flow('vxdata.su','byh src.su c11h c13h c55h data.su', ''' /bin/cp ${SOURCES[5]} ${TARGETS[0]} && ''' + esg + ''' grid_info=${SOURCES[0]} buoyancy=${SOURCES[0]} c11=${SOURCES[2]} c22=${SOURCES[2]} c33=${SOURCES[2]} c23=${SOURCES[3]} c13=${SOURCES[3]} c12=${SOURCES[3]} c44=${SOURCES[4]} c55=${SOURCES[4]} c66=${SOURCES[4]} source_vx=${SOURCES[1]} data_vx=${TARGETS[0]} deriv=0 adjoint=0 order=3 cfl=0.5 cs=%g cp=%g dump_ldc=1 dt=%g eflag=0 '''%(cs,cp,dt), stdin=0, stdout=0) ## (2) create esg movie Flow('movie','byh src.su c11h c13h c55h', ''' spike n1=%d n2=%d n3=%d n4=60 d1=%d d2=%d d3=%d d4=5 label1="Distance" unit1="m" label2="Distance" unit2="m" label3="Distance" unit3="m" label4="Time" unit4="ms" | scale dscale=0 | put dim=3 gdim=4 id1=0 id2=1 id3=2 id4=3 > ${TARGETS[0]} && '''%(nz,nx,ny,dz,dx,dy) + esg + ''' grid_info=${SOURCES[0]} buoyancy=${SOURCES[0]} c11=${SOURCES[2]} c22=${SOURCES[2]} c33=${SOURCES[2]} c23=${SOURCES[3]} c13=${SOURCES[3]} c12=${SOURCES[3]} c44=${SOURCES[4]} c55=${SOURCES[4]} c66=${SOURCES[4]} source_vx=${SOURCES[1]} movie_vx=${TARGETS[0]} deriv=0 adjoint=0 order=3 cfl=0.5 cs=%g cp=%g dump_ldc=1 dt=%g eflag=0 '''%(cs,cp,dt), stdin=0, stdout=0) ## (3) check energy trace Flow('EnergyTrace','byh src.su c11h c13h c55h', ''' spike n1=%d n2=%d n3=%d n4=120 d1=%d d2=%d d3=%d d4=5 label1="Distance" unit1="m" label2="Distance" unit2="m" label3="Distance" unit3="m" label4="Time" unit4="ms" | scale dscale=0 | put dim=3 gdim=4 id1=0 id2=1 id3=2 id4=3 > ${TARGETS[0]} && '''%(nz,nx,ny,dz,dx,dy) + esg + ''' grid_info=${SOURCES[0]} buoyancy=${SOURCES[0]} c11=${SOURCES[2]} c22=${SOURCES[2]} c33=${SOURCES[2]} c23=${SOURCES[3]} c13=${SOURCES[3]} c12=${SOURCES[3]} c44=${SOURCES[4]} c55=${SOURCES[4]} c66=${SOURCES[4]} source_vx=${SOURCES[1]} data_trace=${TARGETS[0]} deriv=0 adjoint=0 order=3 cfl=0.5 cs=%g cp=%g dump_ldc=1 dt=%g eflag=1 '''%(cs,cp,dt), stdin=0, stdout=0) ### === PREPROCESSING # movie.vpl shows the 3D movie Plot('movie', ''' byte gainpanel=all | grey4 frame1=%d frame2=%d frame3=%d clip=0.01 wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n '''%(nz/2,nx/2,ny/2)) # trace.vpl shows the e(t) Flow('trace','EnergyTrace','window n1=1 n2=1 n3=1 f1=1 f2=1 f3=1') Plot('trace','sfgraph') End()