Skip to content

Commit b85b2e3

Browse files
author
George Bisbas
committed
examples: Update 13/15/16
1 parent 394994b commit b85b2e3

3 files changed

Lines changed: 69 additions & 73 deletions

File tree

‎examples/seismic/tutorials/13_LSRTM_acoustic.ipynb‎

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -469,41 +469,40 @@
469469
"image = np.zeros((model0.vp.shape[0], model0.vp.shape[1]))\n",
470470
"\n",
471471
"nrec=101\n",
472-
"niter=20 # number of iterations of the LSRTM\n",
473-
"history = np.zeros((niter, 1)) #objective function\n",
472+
"niter=20 # Number of iterations of the LSRTM\n",
473+
"history = np.zeros((niter, 1)) # Objective function\n",
474474
"\n",
475475
"image_prev = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
476-
" \n",
476+
"\n",
477477
"grad_prev = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
478478
"\n",
479479
"yk = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
480480
"\n",
481481
"sk = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
482482
"\n",
483-
"for k in range(niter) :\n",
484-
" \n",
485-
" dm = image_up_dev # Reflectivity for Calculated data via Born\n",
483+
"for k in range(niter):\n",
484+
"\n",
485+
" dm = image_up_dev # Reflectivity for Calculated data via Born\n",
486486
"\n",
487487
" print('LSRTM Iteration',k+1)\n",
488-
" \n",
488+
"\n",
489489
" objective,grad_full,d_obs,d_syn = lsrtm_gradient(dm)\n",
490-
" \n",
490+
"\n",
491491
" history[k] = objective\n",
492-
" \n",
492+
"\n",
493493
" yk = grad_full.data - grad_prev\n",
494-
" \n",
494+
"\n",
495495
" sk = image_up_dev - image_prev\n",
496496
"\n",
497497
" alfa = get_alfa(yk,sk,k)\n",
498-
" \n",
498+
"\n",
499499
" grad_prev = grad_full.data\n",
500500
"\n",
501501
" image_prev = image_up_dev\n",
502-
" \n",
502+
"\n",
503503
" image_up_dev = image_up_dev - alfa*grad_full.data\n",
504-
" \n",
504+
"\n",
505505
" if k == 0: # Saving the first migration using Born operator.\n",
506-
" \n",
507506
" image = image_up_dev"
508507
]
509508
},
@@ -560,7 +559,7 @@
560559
" vmin=-.05,\n",
561560
" vmax=.05,\n",
562561
" cmap=cmap,extent=extent)\n",
563-
" \n",
562+
"\n",
564563
" plt.xlabel('X position (km)')\n",
565564
" plt.ylabel('Depth (km)')\n",
566565
"\n",
@@ -570,6 +569,7 @@
570569
" divider = make_axes_locatable(ax)\n",
571570
" cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
572571
" plt.colorbar(plot, cax=cax)\n",
572+
"\n",
573573
" plt.show()\n"
574574
]
575575
},
@@ -600,7 +600,7 @@
600600
],
601601
"source": [
602602
"#NBVAL_IGNORE_OUTPUT\n",
603-
"slices=tuple(slice(model.nbl,-model.nbl) for _ in range(2))\n",
603+
"slices=tuple(slice(model.nbl, -model.nbl) for _ in range(2))\n",
604604
"rtm = image[slices]\n",
605605
"plot_image(np.diff(rtm, axis=1))"
606606
]

‎examples/seismic/tutorials/15_tti_qp_pure.ipynb‎

Lines changed: 30 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@
4747
"source": [
4848
"import numpy as np\n",
4949
"from devito import (Function, TimeFunction, cos, sin, solve,\n",
50-
" Eq, Operator, configuration, norm)\n",
50+
" Eq, Operator)\n",
5151
"from examples.seismic import TimeAxis, RickerSource, Receiver, demo_model\n",
5252
"from matplotlib import pyplot as plt"
5353
]
@@ -65,34 +65,23 @@
6565
"execution_count": 2,
6666
"id": "4f545ff1",
6767
"metadata": {},
68-
"outputs": [
69-
{
70-
"name": "stderr",
71-
"output_type": "stream",
72-
"text": [
73-
"Operator `pad_vp` ran in 0.01 s\n",
74-
"Operator `pad_epsilon` ran in 0.01 s\n",
75-
"Operator `pad_delta` ran in 0.01 s\n",
76-
"Operator `pad_theta` ran in 0.01 s\n"
77-
]
78-
}
79-
],
68+
"outputs": [],
8069
"source": [
81-
"# NBVAL_IGNORE_OUTPUT \n",
70+
"# NBVAL_IGNORE_OUTPUT\n",
8271
"\n",
83-
"shape = (101,101) # 101x101 grid\n",
84-
"spacing = (10.,10.) # spacing of 10 meters\n",
85-
"origin = (0.,0.) \n",
72+
"shape = (101,101) # 101x101 grid\n",
73+
"spacing = (10.,10.) # spacing of 10 meters\n",
74+
"origin = (0.,0.)\n",
8675
"nbl = 0 # number of pad points\n",
8776
"\n",
8877
"model = demo_model('layers-tti', spacing=spacing, space_order=8,\n",
8978
" shape=shape, nbl=nbl, nlayers=1)\n",
9079
"\n",
9180
"# initialize Thomsem parameters to those used in Mu et al., (2020)\n",
92-
"model.update('vp', np.ones(shape)*3.6) # km/s\n",
81+
"model.update('vp', np.ones(shape)*3.6) # km/s\n",
9382
"model.update('epsilon', np.ones(shape)*0.23)\n",
9483
"model.update('delta', np.ones(shape)*0.17)\n",
95-
"model.update('theta', np.ones(shape)*(45.*(np.pi/180.))) # radians"
84+
"model.update('theta', np.ones(shape)*(45.*(np.pi/180.))) # radians"
9685
]
9786
},
9887
{
@@ -233,7 +222,7 @@
233222
"bc += [Eq(pp[t+1,0, z], 0.)]\n",
234223
"bc += [Eq(pp[t+1,shape[0]-1+2*nbl, z], 0.)]\n",
235224
"\n",
236-
"# set source and receivers\n",
225+
"# Set source and receivers\n",
237226
"src = RickerSource(name='src',grid=model.grid,f0=0.02,npoint=1,time_range=time_range)\n",
238227
"src.coordinates.data[:,0] = model.domain_size[0]* .5\n",
239228
"src.coordinates.data[:,1] = model.domain_size[0]* .5\n",
@@ -243,14 +232,18 @@
243232
"rec = Receiver(name='rec',grid=model.grid,npoint=shape[0],time_range=time_range)\n",
244233
"rec.coordinates.data[:, 0] = np.linspace(model.origin[0],model.domain_size[0], num=model.shape[0])\n",
245234
"rec.coordinates.data[:, 1] = 2*spacing[1]\n",
235+
"\n",
246236
"# Create interpolation expression for receivers\n",
247237
"rec_term = rec.interpolate(expr=p.forward)\n",
248238
"\n",
249239
"# Operators\n",
250-
"optime=Operator([update_p] + src_term + rec_term)\n",
251-
"oppres=Operator([update_q] + bc)\n",
240+
"optime = Operator([update_p] + src_term + rec_term)\n",
241+
"oppres = Operator([update_q] + bc)\n",
242+
"\n",
252243
"\n",
253-
"# you can print the generated code for both operators by typing print(optime) and print(oppres)"
244+
"# You can print the generated code for both operators by uncommenting the following lines\n",
245+
"# print(optime)\n",
246+
"# print(oppres)"
254247
]
255248
},
256249
{
@@ -620,17 +613,17 @@
620613
],
621614
"source": [
622615
"# NBVAL_IGNORE_OUTPUT\n",
623-
"psave =np.empty ((time_range.num,model.grid.shape[0],model.grid.shape[1]))\n",
616+
"psave =np.empty((time_range.num,model.grid.shape[0], model.grid.shape[1]))\n",
624617
"niter_poisson = 1200\n",
625618
"\n",
626619
"# This is the time loop.\n",
627-
"for step in range(0,time_range.num-2):\n",
628-
" q.data[:,:]=pp.data[(niter_poisson+1)%2,:,:]\n",
620+
"for step in range(0, time_range.num-2):\n",
621+
" q.data[:, :]=pp.data[(niter_poisson+1)%2, :, :]\n",
629622
" optime(time_m=step, time_M=step, dt=dt)\n",
630-
" pp.data[:,:]=0.\n",
631-
" b.data[:,:]=p.data[(step+1)%3,:,:]\n",
623+
" pp.data[:, :]=0.\n",
624+
" b.data[:, :]=p.data[(step+1)%3, :, :]\n",
632625
" oppres(time_M = niter_poisson)\n",
633-
" psave[step,:,:]=p.data[(step+1)%3,:,:]"
626+
" psave[step, :, :]=p.data[(step+1)%3, :, :]"
634627
]
635628
},
636629
{
@@ -691,18 +684,18 @@
691684
"fig, axes = plt.subplots(2, 5, figsize=(18, 7), sharex=True)\n",
692685
"fig.suptitle(\"Snapshots\", size=14)\n",
693686
"for count, ax in enumerate(axes.ravel()):\n",
694-
" snapshot = factor*count\n",
687+
" snapshot = factor * count\n",
695688
" ax.imshow(np.transpose(psave[snapshot,:,:]), cmap=\"seismic\",\n",
696-
" vmin=-amax, vmax=+amax, extent=plt_extent)\n",
697-
" ax.plot(model.domain_size[0]* .5, model.domain_size[1]* .5, \\\n",
698-
" 'red', linestyle='None', marker='*', markersize=8, label=\"Source\")\n",
689+
" vmin=-amax, vmax=+amax, extent=plt_extent)\n",
690+
" ax.plot(model.domain_size[0] * 0.5, model.domain_size[1] * 0.5,\n",
691+
" 'red', linestyle='None', marker='*', markersize=8, label=\"Source\")\n",
699692
" ax.grid()\n",
700-
" ax.tick_params('both', length=2, width=0.5, which='major',labelsize=10)\n",
701-
" ax.set_title(\"Wavefield at t=%.2fms\" % (factor*count*dt),fontsize=10)\n",
693+
" ax.tick_params('both', length=2, width=0.5, which='major', labelsize=10)\n",
694+
" ax.set_title(f\"Wavefield at t={factor * count * dt:.2f}ms\", fontsize=10)\n",
702695
"for ax in axes[1, :]:\n",
703-
" ax.set_xlabel(\"X Coordinate (m)\",fontsize=10)\n",
696+
" ax.set_xlabel(\"X Coordinate (m)\", fontsize=10)\n",
704697
"for ax in axes[:, 0]:\n",
705-
" ax.set_ylabel(\"Z Coordinate (m)\",fontsize=10)"
698+
" ax.set_ylabel(\"Z Coordinate (m)\", fontsize=10)"
706699
]
707700
},
708701
{

‎examples/seismic/tutorials/16_ader_fd.ipynb‎

Lines changed: 23 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -49,11 +49,12 @@
4949
"source": [
5050
"#NBVAL_IGNORE_OUTPUT\n",
5151
"# Necessary imports\n",
52-
"import devito as dv\n",
5352
"import sympy as sp\n",
5453
"import numpy as np\n",
5554
"import matplotlib.pyplot as plt\n",
5655
"\n",
56+
"from devito import (Function, TimeFunction, Grid, Operator, Eq, solve,\n",
57+
" VectorTimeFunction, grad, div)\n",
5758
"from examples.seismic import TimeAxis, RickerSource"
5859
]
5960
},
@@ -73,9 +74,9 @@
7374
"outputs": [],
7475
"source": [
7576
"# 1km x 1km grid\n",
76-
"grid = dv.Grid(shape=(201, 201), extent=(1000., 1000.))\n",
77-
"p = dv.TimeFunction(name='p', grid=grid, space_order=16)\n",
78-
"v = dv.VectorTimeFunction(name='v', grid=grid, space_order=16, staggered=(None, None))"
77+
"grid = Grid(shape=(201, 201), extent=(1000., 1000.))\n",
78+
"p = TimeFunction(name='p', grid=grid, space_order=16)\n",
79+
"v = VectorTimeFunction(name='v', grid=grid, space_order=16, staggered=(None, None))"
7980
]
8081
},
8182
{
@@ -94,8 +95,8 @@
9495
"outputs": [],
9596
"source": [
9697
"# Material parameters\n",
97-
"c = dv.Function(name='c', grid=grid)\n",
98-
"rho = dv.Function(name='rho', grid=grid)\n",
98+
"c = Function(name='c', grid=grid)\n",
99+
"rho = Function(name='rho', grid=grid)\n",
99100
"\n",
100101
"c.data[:] = 1.5\n",
101102
"c.data[:, :150] = 1.25\n",
@@ -126,8 +127,8 @@
126127
"metadata": {},
127128
"outputs": [],
128129
"source": [
129-
"# dv.grad(dv.div(v)) is not the same as expanding the continuous operator and then discretising\n",
130-
"# This is because dv.grad(dv.div(v)) applies a gradient stencil to a divergence stencil\n",
130+
"# grad(div(v)) is not the same as expanding the continuous operator and then discretising\n",
131+
"# This is because grad(div(v)) applies a gradient stencil to a divergence stencil\n",
131132
"def graddiv(f):\n",
132133
" return sp.Matrix([[f[0].dx2 + f[1].dxdy],\n",
133134
" [f[0].dxdy + f[1].dy2]])\n",
@@ -147,8 +148,8 @@
147148
" return f.dx4 + 2*f.dx2dy2 + f.dy4\n",
148149
"\n",
149150
"# First time derivatives\n",
150-
"pdt = rho*c2*dv.div(v)\n",
151-
"vdt = b*dv.grad(p)\n",
151+
"pdt = rho*c2*div(v)\n",
152+
"vdt = b*grad(p)\n",
152153
"\n",
153154
"# Second time derivatives\n",
154155
"pdt2 = c2*p.laplace\n",
@@ -200,8 +201,8 @@
200201
"dt = grid.stepping_dim.spacing\n",
201202
"\n",
202203
"# Update equations (4th-order ADER timestepping)\n",
203-
"eq_p = dv.Eq(p.forward, p + dt*pdt + (dt**2/2)*pdt2 + (dt**3/6)*pdt3 + (dt**4/24)*pdt4)\n",
204-
"eq_v = dv.Eq(v.forward, v + dt*vdt + (dt**2/2)*vdt2 + (dt**3/6)*vdt3 + (dt**4/24)*vdt4)"
204+
"eq_p = Eq(p.forward, p + dt*pdt + (dt**2/2)*pdt2 + (dt**3/6)*pdt3 + (dt**4/24)*pdt4)\n",
205+
"eq_v = Eq(v.forward, v + dt*vdt + (dt**2/2)*vdt2 + (dt**3/6)*vdt3 + (dt**4/24)*vdt4)"
205206
]
206207
},
207208
{
@@ -286,7 +287,7 @@
286287
"#NBVAL_IGNORE_OUTPUT\n",
287288
"src_term = src.inject(field=p.forward, expr=src)\n",
288289
"\n",
289-
"op = dv.Operator([eq_p, eq_v] + src_term)\n",
290+
"op = Operator([eq_p, eq_v] + src_term)\n",
290291
"op.apply(dt=op_dt)"
291292
]
292293
},
@@ -355,20 +356,22 @@
355356
}
356357
],
357358
"source": [
359+
"from devito.types import NODE\n",
360+
"\n",
358361
"#NBVAL_IGNORE_OUTPUT\n",
359-
"ps = dv.TimeFunction(name='ps', grid=grid, space_order=16, staggered=dv.NODE)\n",
360-
"vs = dv.VectorTimeFunction(name='vs', grid=grid, space_order=16)\n",
362+
"ps = TimeFunction(name='ps', grid=grid, space_order=16, staggered=NODE)\n",
363+
"vs = VectorTimeFunction(name='vs', grid=grid, space_order=16)\n",
361364
"\n",
362365
"# First time derivatives\n",
363-
"psdt = rho*c2*dv.div(vs.forward)\n",
364-
"vsdt = b*dv.grad(ps)\n",
366+
"psdt = rho*c2*div(vs.forward)\n",
367+
"vsdt = b*grad(ps)\n",
365368
"\n",
366-
"eq_ps = dv.Eq(ps.forward, ps + dt*psdt)\n",
367-
"eq_vs = dv.Eq(vs.forward, vs + dt*vsdt)\n",
369+
"eq_ps = Eq(ps.forward, ps + dt*psdt)\n",
370+
"eq_vs = Eq(vs.forward, vs + dt*vsdt)\n",
368371
"\n",
369372
"src_term_s = src.inject(field=ps.forward, expr=src)\n",
370373
"\n",
371-
"ops = dv.Operator([eq_vs, eq_ps] + src_term_s)\n",
374+
"ops = Operator([eq_vs, eq_ps] + src_term_s)\n",
372375
"ops.apply(dt=op_dt)"
373376
]
374377
},

0 commit comments

Comments
 (0)