from simpegseis import *
%pylab inline
from SimPEG.Utils import mkvc
Populating the interactive namespace from numpy and matplotlib
# Step1: set time
time = np.linspace(0, 0.04, 2**9)
dt = time[1]-time[0]
# Step2: set Tx and Rx
options={'tlag':0.0025, 'fmain':400} # You need to set waveform to set Tx
rx = AcousticRx(np.vstack((np.r_[0, 1], np.r_[0, 1])))
tx = AcousticTx(np.r_[0, 1], time, [rx], **options)
# Step3: set survey (pass txlist)
survey = SurveyAcoustic([tx])
wave = tx.RickerWavelet()
# Step4: set mesh
cs = 0.5
hx = np.ones(150)*cs
hy = np.ones(150)*cs
mesh = Mesh.TensorMesh([hx, hy], 'CC')
# Step5: set problem (pass mesh) and pair with survey
prob = AcousticProblemPML(mesh)
prob.pair(survey)
# Step6: set boundary
prob.setPMLBC(20, dt, bcflag='all', const=0.)
prob.storefield = True
# Step7: set velocity model and check stability
v = np.ones(mesh.nC)*2000.
prob.stabilitycheck(v, time, 100.)
You are good to go:) >> Stability information dt: 7.83e-05 s Optimal dt: 1.25e-04 s Cell per wavelength (G): 4.00e+01 Optimal G: 1.60e+01
# Step8: run forward
U = prob.fields(v)
>> Start Computing Acoustic Wave >> dt: 7.83e-05 s >> Optimal dt: 1.25e-04 s >> Main frequency, fmain: 1.00e+02 Hz >> Cell per wavelength (G): 4.00e+01 Tx at ( 0.00, 0.00): 1/ 1 >>Elapsed time: 2.55e+00 s
# Step9: project data
data = survey.projectFields(U)
print data[0].shape
(2L, 512L)
from JSAnimation import IPython_display
from matplotlib import animation
fig, ax = plt.subplots(1,1, figsize = (8, 8))
ax.set_xlabel('Easting (m)', fontsize = 16)
ax.set_ylabel('Depth (m)', fontsize = 16)
extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]
ax.set_xlim(extent[:2])
ax.set_ylim(extent[2:])
nskip = 20
def animate(i_id):
icount = i_id*nskip
frame = ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent)
return frame
animation.FuncAnimation(fig, animate, frames=14, interval=40, blit=True)
prob.setPMLBC(30, dt, bcflag='all', const=1.)
U = prob.fields(v)
>> Start Computing Acoustic Wave >> dt: 7.83e-05 s >> Optimal dt: 1.25e-04 s >> Main frequency, fmain: 1.00e+02 Hz >> Cell per wavelength (G): 4.00e+01 Tx at ( 0.00, 0.00): 1/ 1 >>Elapsed time: 2.75e+00 s
from JSAnimation import IPython_display
from matplotlib import animation
fig, ax = plt.subplots(1,1, figsize = (8, 8))
ax.set_xlabel('Easting (m)', fontsize = 16)
ax.set_ylabel('Depth (m)', fontsize = 16)
extent = [mesh.vectorCCx.min(), mesh.vectorCCx.max(), mesh.vectorCCy.min()-mesh.vectorCCy.max(), mesh.vectorCCy.max()-mesh.vectorCCy.max()]
ax.set_xlim(extent[:2])
ax.set_ylim(extent[2:])
nskip = 20
def animate(i_id):
icount = i_id*nskip
frame = ax.imshow(np.flipud(U[0][:,icount].reshape((mesh.nCx, mesh.nCy), order = 'F').T), cmap = 'binary', extent=extent)
return frame
animation.FuncAnimation(fig, animate, frames=14, interval=40, blit=True)
# !ipython nbconvert --to html SeimsicFWD.ipynb