#!/bin/python # Sconstruct for 2D 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: # Somehow if you are running this script on MACBOOK, you have to run "scons" # twice. # If you want to run it on Stampede, make sure the email account in "penv" # is configured to your own account, and "SetEnv" is set correctly. Shown below # provides you two "SetEnv"s for Stampede and DAVinCI respectively. #SetEnv = {'platf':'stampede','nodes':str(1),'ppn':str(NP1*NP2),'wall':'00:20:00'} #SetEnv = {'platf':'mpi', 'ppn':str(NP1*NP2)} from rsf.proj import * from newbatch import tripExec import os from math import ceil 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/svn-esg/main/esg.x') fetches = {} ### === SET PARAMETERS import math ## for source fpeak=20 T=0.2 # source duration ## for stiffness fields rho=1 cs=1 cp=1.5 mu0=rho*cs**2; lambda0=rho*cp**2-2*mu0; ## for grid nz=401 nx=nz-30 dz=cs*1000/fpeak/30 #unit[m] dx=0.9*dz dt=min(dz,dx)/cp/math.sqrt(2)*0.5 #unit[ms] ### === CREATE FIELDS Flow('byh', None, ''' makevel n1=%d n2=%d d1=%f d2=%f v000=1 | sfput dim=2 gdim=2 id1=0 id2=1 '''%(nz,nx,dz,dx), stdin=0) Flow('c11h-1', None, ''' makevel n1=%d n2=%d d1=%f d2=%f v000=%f | sfput dim=2 gdim=2 id1=0 id2=1 '''%(ceil(nz*0.7),nx,dz,dx, lambda0), stdin=0) Flow('c11h-2', None, ''' makevel n1=%d n2=%d d1=%f d2=%f v000=%f | sfput dim=2 gdim=2 id1=0 id2=1 '''%(nz-ceil(nz*0.7),nx,dz,dx, lambda0+2*mu0), stdin=0) Flow('c11h',['c11h-1', 'c11h-2'],'''cat ${SOURCES[1:-1]} axis=1''') Flow('c13h', None, ''' makevel n1=%d n2=%d d1=%f d2=%f v000=%f | sfput dim=2 gdim=2 id1=0 id2=1 '''%(nz,nx,dz,dx, lambda0), stdin=0) Flow('c55h-1', None, ''' makevel n1=%d n2=%d d1=%f d2=%f v000=%f | sfput dim=2 gdim=2 id1=0 id2=1 '''%(ceil(nz*0.7),nx,dz,dx, 0), stdin=0) Flow('c55h-2', None, ''' makevel n1=%d n2=%d d1=%f d2=%f v000=%f | sfput dim=2 gdim=2 id1=0 id2=1 '''%(nz-ceil(nz*0.7),nx,dz,dx, mu0), stdin=0) Flow('c55h',['c55h-1', 'c55h-2'],'''cat ${SOURCES[1:-1]} axis=1''') ### === CREATE SOURCE (time unit[s]) 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=%d | '%(ceil(dx*nx/2)) + sushw + ' key=gelev a=%d '%(-ceil(dz*nz/2)), stdin=0) Flow('src_data.su',None, sunull + ' nt=300 ntr=1 | ' + sushw + ' key=sx a=%d b=%d | '%(0, 100) + sushw + ' key=selev a=0', stdin=0) Flow('src_towed.su', ['src.su', 'src_data.su'], '''towed_array src=${SOURCES[0]} data=${SOURCES[1]} towed=${TARGETS[0]}''', stdin=0, stdout=0) ### === CREATE RECEIVER TRACE FILE (time unit[s]) Flow('data.su', None, sunull + ' nt=1001 ntr=%d dt=0.002 | '%(ceil(nx/2)) + sushw + ' key=gx a=%d b=%d | '%(ceil(nx/2*dx)-ceil(dx)*ceil(nx/4),ceil(dx)) + sushw + ' key=gelev a=%d > ${TARGETS[0]}'%(-ceil(dz*nz/3)), stdin=0) ### === CREATE MOVIE FILE (time unit[ms]) Flow('movie', None, '''makevel n1=%d n2=%d n3=150 d1=%f d2=%f d3=10 v000=0.0 | sfput dim=2 gdim=3 id1=0 id2=1 id3=2 '''%(nz,nx,dz,dx), stdin=0) ### === START jobs=[] # create vx(/vz)data, vx(/vz)movie profile for non-effective media(0) # and effective(1) media for i in range(2): for NP1 in [2]: for NP2 in [2]: jobs = jobs + [{'job': 'esg-data' + str(NP1) + str(NP2), 'pre': ''' /bin/cp ${SOURCES[5]} ${TARGETS[0]} && /bin/cp ${SOURCES[5]} ${TARGETS[1]} && cp ${SOURCES[6]} ${TARGETS[2]} && cp ${SOURCES[6]} ${TARGETS[3]} ''', 'cmd': esg + ''' grid_info=${SOURCES[0]} buoyancy=${SOURCES[0]} c11=${SOURCES[1]} c33=${SOURCES[1]} c13=${SOURCES[2]} c55=${SOURCES[3]} source_sxx=${SOURCES[4]} source_szz=${SOURCES[4]} data_vx=${TARGETS[0]} data_vz=${TARGETS[1]} movie_vx=${TARGETS[2]} movie_vz=${TARGETS[3]} deriv=0 adjoint=0 order=2 dump_ldc=1 dt=%g eflag=0 mpi_np1=%d mpi_np2=%d effective_media=%d '''%(dt,NP1,NP2,i), 'src':['byh.rsf', 'c11h.rsf', 'c13h.rsf', 'c55h.rsf', 'src_towed.su', 'data.su', 'movie'], 'tgt':['vxdata'+str(i)+'.su', 'vzdata'+str(i)+'.su', 'vxmovie'+str(i)+'.rsf', 'vzmovie'+str(i)+'.rsf'], 'exe':{'platf':'mpi', 'ppn':str(NP1*NP2)} },] Flow(['vxdata'+str(i), 'tfile-x'+str(i)],'vxdata'+str(i)+'.su','''suread tfile=${TARGETS[1]} endian=0''') Flow(['vzdata'+str(i), 'tfile-z'+str(i)],'vzdata'+str(i)+'.su','''suread tfile=${TARGETS[1]} endian=0''') Plot('vxdata'+str(i), '''sfgrey clip=0.01 scalebar=y color=e title=\"vx%d\" '''%(i)) Plot('vzdata'+str(i), '''sfgrey clip=0.01 scalebar=y color=e title=\"vz%d\" '''%(i)) Plot('vxmovie'+str(i), 'window f3=61 n3=1 | sfgrey clip=0.01 color=e scalebar=y') Plot('vzmovie'+str(i), 'window f3=61 n3=1 | sfgrey clip=0.01 color=e scalebar=y') tripExec(jobs,penv) ### === PREPROCESSING (COMPARE DIFFERENCE) Flow('vxdata-diff',['vxdata0', 'vxdata1'], '''math x=${SOURCES[0]} y=${SOURCES[1]} output='x-y' ''',stdin=0) Plot('vxdata-diff','''grey clip=0.01 scalebar=y color=e title=\"vx-diff\" ''') Flow('vzdata-diff', ['vzdata0', 'vzdata1'], '''math x=${SOURCES[0]} y=${SOURCES[1]} output='x-y' ''',stdin=0) Plot('vzdata-diff','''grey clip=0.01 scalebar=y color=e title=\"vz-diff\" ''') Flow('vxmovie-diff', ['vxmovie0', 'vxmovie1'], '''math x=${SOURCES[0]} y=${SOURCES[1]} output='x-y' ''',stdin=0) Plot('vxmovie-diff','''window f3=61 n3=1 | grey clip=0.01 scalebar=y color=e title=\"vxmovie-diff\" ''') Flow('vzmovie-diff', ['vzmovie0', 'vzmovie1'], '''math x=${SOURCES[0]} y=${SOURCES[1]} output='x-y' ''',stdin=0) Plot('vzmovie-diff','''window f3=61 n3=1 | grey clip=0.01 scalebar=y color=e title=\"vzmovie-diff\" ''') Flow('vxmovie-add', ['vxmovie0', 'vxmovie1'], '''math x=${SOURCES[0]} y=${SOURCES[1]} output='x+y' ''',stdin=0) Plot('vxmovie-add','''window f3=61 n3=1 | grey clip=0.01 scalebar=y color=e title=\"vxmovie-add\" ''') Flow('vzmovie-add', ['vzmovie0', 'vzmovie1'], '''math x=${SOURCES[0]} y=${SOURCES[1]} output='x+y' ''',stdin=0) Plot('vzmovie-add','''window f3=61 n3=1 | grey clip=0.01 scalebar=y color=e title=\"vzmovie-add\" ''') End()