# Sconstruct for 3D esg: (1) record trace data; # (2) generate movie; # (3) display discrete total energy of the solution as a function of time # (parallel version) # Author: MZ # # Note: # (1) Somehow on Macbook the following three jobs in this script cannot be executed by doing scons # once, you need to run scons three times in order to execute 3 jobs. You should see errors after # running scons for the first and second time. But don't be panic, keep running scons will yield # desired results. If some warmhearted person figure out how to fix this problem -- run scons once and # run three jobs, send me an email (mz10@rice.edu). I will be very much appreciated. # (2) If you are running this script on computing platforms (e.g. stampede, Davinci), please change # the email setting in the following penv. # (3) Currently this script runs on Mac, if you want to run it on Stampede, # copy one of the following lines to jobs:'exe'. #SetEnv = {'platf':'stampede','nodes':str(1),'ppn':str(NP1*NP2*NP3),'wall':'00:10:00'} #SetEnv = {'platf':'mpi', 'ppn':str(NP1*NP2*NP3)} from rsf.proj import * from newbatch import tripExec import os penv = {'stampede' : { 'batch' : 'slurm', 'queue' : 'normal', 'acct' : 'FDTD3D-Cont', 'mail' : '***@rice.edu', 'launcher' : 'ibrun'}, 'davinci' : { 'batch' : 'slurm', 'queue' : 'trip', 'acct' : '***', 'mail' : '***@rice.edu', 'launcher' : 'ibrun'} } 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') MYAPPS = os.getenv('RSFSRC') towed_array = os.path.join(MYAPPS,'trip2.1/iwave/trace/main/towed_array.x') stdmdl = os.path.join(MYAPPS,'trip2.1/iwave/grid/main/standardmodel.x') esg = os.path.join(MYAPPS,'trip2.1/iwave/esg/main/esg.x') fetches = {} ### === 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 jobs=[] ## (1) create data_vx profile for NP1 in [1]: for NP2 in [2]: for NP3 in [2]: jobs = jobs + [{'job': 'esg-data' + str(NP1) + str(NP2) + str(NP3), 'pre': ''' /bin/cp ${SOURCES[5]} ${TARGETS[0]} ''', 'cmd': 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.4 cs=%g cp=%g dt=%g dump_ldc=1 eflag=0 mpi_np1=%d mpi_np2=%d mpi_np3=%d '''%(cs,cp,dt,NP1,NP2,NP3), 'src':['byh.rsf', 'src.su', 'c11h.rsf', 'c13h.rsf', 'c55h.rsf','data.su'], 'tgt':['vxdata.su'], 'exe':{'platf':'mpi', 'ppn':str(NP1*NP2*NP3)} },] ## (2) create esg movie for NP1 in [1]: for NP2 in [2]: for NP3 in [2]: jobs = jobs + [{'job': 'esg-sim' + str(NP1) + str(NP2) + str(NP3), 'pre': ''' 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), 'cmd': 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.4 cs=%g cp=%g dump_ldc=1 dt=%g dump_ldc=1 eflag=0 mpi_np1=%d mpi_np2=%d mpi_np3=%d '''%(cs,cp,dt,NP1,NP2,NP3), 'src':['byh.rsf', 'src.su', 'c11h.rsf', 'c13h.rsf', 'c55h.rsf'], 'tgt':['movie.rsf'], 'exe':{'platf':'mpi', 'ppn':str(NP1*NP2*NP3)} },] ## (3) check energy trace for NP1 in [1]: for NP2 in [2]: for NP3 in [2]: jobs = jobs + [{'job': 'esg-energy' + str(NP1) + str(NP2) + str(NP3), 'pre': ''' 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), 'cmd': 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.4 cs=%g cp=%g dt=%g dump_ldc=1 eflag=1 mpi_np1=%d mpi_np2=%d mpi_np3=%d '''%(cs,cp,dt,NP1,NP2,NP3), 'src':['byh.rsf', 'src.su', 'c11h.rsf', 'c13h.rsf', 'c55h.rsf'], 'tgt':['EnergyTrace.rsf'], 'exe':{'platf':'mpi', 'ppn':str(NP1*NP2*NP3)} },] tripExec(jobs,penv) ### === 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()