Skip to content

Commit 5fec3ad

Browse files
George Bisbasgeorgebisbas
George Bisbas
authored and
georgebisbas
committed
examples: Further reviews and updates
1 parent e4736a9 commit 5fec3ad

2 files changed

Lines changed: 50 additions & 58 deletions

File tree

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

Lines changed: 45 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -126,27 +126,18 @@
126126
"%matplotlib inline\n",
127127
"import numpy as np\n",
128128
"\n",
129-
"from devito import Operator,Eq,solve,Grid,SparseFunction,norm\n",
130-
"from devito import TimeFunction,Function\n",
131-
"from devito import gaussian_smooth\n",
132-
"from devito import mmax\n",
129+
"from devito import (Operator, Eq, solve, Grid, norm, Function,\n",
130+
" gaussian_smooth, mmax)\n",
133131
"\n",
134-
"from devito.logger import info\n",
135-
"\n",
136-
"from examples.seismic import Model\n",
137-
"from examples.seismic import plot_velocity,plot_shotrecord\n",
138-
"from examples.seismic import Receiver\n",
139-
"from examples.seismic import PointSource\n",
140-
"from examples.seismic import plot_image,AcquisitionGeometry\n",
141-
"from examples.seismic import TimeAxis\n",
132+
"from examples.seismic import (Model, plot_velocity, Receiver, plot_image,\n",
133+
" AcquisitionGeometry, TimeAxis)\n",
142134
"\n",
143135
"from examples.seismic.self_adjoint import (setup_w_over_q)\n",
144136
"from examples.seismic.acoustic import AcousticWaveSolver\n",
145137
"\n",
146138
"import matplotlib.pyplot as plt\n",
147139
"from mpl_toolkits.axes_grid1 import ImageGrid\n",
148140
"from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable\n",
149-
"import matplotlib.ticker as plticker\n",
150141
"\n",
151142
"from devito import configuration\n",
152143
"configuration['log-level'] = 'WARNING'"
@@ -183,7 +174,7 @@
183174
"omega = 2.0 * np.pi * fpeak\n",
184175
"qmin = 0.1\n",
185176
"qmax = 100000\n",
186-
"npad=50\n",
177+
"npad = 50\n",
187178
"dtype = np.float32\n",
188179
"\n",
189180
"nshots = 21\n",
@@ -309,57 +300,57 @@
309300
"outputs": [],
310301
"source": [
311302
"def lsrtm_gradient(dm):\n",
312-
" \n",
303+
"\n",
313304
" residual = Receiver(name='residual', grid=model.grid, time_range=geometry.time_axis,\n",
314305
" coordinates=geometry.rec_positions)\n",
315-
" \n",
306+
"\n",
316307
" d_obs = Receiver(name='d_obs', grid=model.grid,time_range=geometry.time_axis,\n",
317308
" coordinates=geometry.rec_positions)\n",
318309
"\n",
319310
" d_syn = Receiver(name='d_syn', grid=model.grid,time_range=geometry.time_axis,\n",
320311
" coordinates=geometry.rec_positions)\n",
321-
" \n",
312+
"\n",
322313
" grad_full = Function(name='grad_full', grid=model.grid)\n",
323-
" \n",
314+
"\n",
324315
" grad_illum = Function(name='grad_illum', grid=model.grid)\n",
325-
" \n",
316+
"\n",
326317
" src_illum = Function (name =\"src_illum\", grid = model.grid)\n",
327318
"\n",
328319
" # Using devito's reference of virtual source\n",
329320
" dm_true = (solver.model.vp.data**(-2) - model0.vp.data**(-2))\n",
330-
" \n",
321+
"\n",
331322
" objective = 0.\n",
332323
" u0 = None\n",
333324
" for i in range(nshots):\n",
334-
" \n",
325+
"\n",
335326
" #Observed Data using Born's operator\n",
336327
" geometry.src_positions[0, :] = source_locations[i, :]\n",
337328
"\n",
338329
" _, u0, _ = solver.forward(vp=model0.vp, save=True, u=u0)\n",
339330
" \n",
340331
" _, _, _,_ = solver.jacobian(dm_true, vp=model0.vp, rec = d_obs)\n",
341-
" \n",
332+
"\n",
342333
" #Calculated Data using Born's operator\n",
343334
" solver.jacobian(dm, vp=model0.vp, rec = d_syn)\n",
344-
" \n",
335+
"\n",
345336
" residual.data[:] = d_syn.data[:]- d_obs.data[:]\n",
346-
" \n",
337+
"\n",
347338
" grad_shot,_ = solver.gradient(rec=residual, u=u0, vp=model0.vp)\n",
348-
" \n",
339+
"\n",
349340
" src_illum_upd = Eq(src_illum, src_illum + u0**2)\n",
350341
" op_src = Operator([src_illum_upd])\n",
351342
" op_src.apply()\n",
352-
" \n",
343+
"\n",
353344
" grad_sum = Eq(grad_full, grad_full + grad_shot)\n",
354345
" op_grad = Operator([grad_sum])\n",
355346
" op_grad.apply()\n",
356-
" \n",
347+
"\n",
357348
" objective += .5*norm(residual)**2\n",
358-
" \n",
349+
"\n",
359350
" grad_f = Eq(grad_illum, grad_full/(src_illum+10**-9))\n",
360351
" op_gradf = Operator([grad_f])\n",
361352
" op_gradf.apply()\n",
362-
" \n",
353+
"\n",
363354
" return objective,grad_illum,d_obs,d_syn"
364355
]
365356
},
@@ -394,31 +385,31 @@
394385
"outputs": [],
395386
"source": [
396387
"def get_alfa(grad_iter,image_iter,niter_lsrtm):\n",
397-
" \n",
398-
" \n",
388+
"\n",
389+
"\n",
399390
" term1 = np.dot(image_iter.reshape(-1), image_iter.reshape(-1))\n",
400-
" \n",
391+
"\n",
401392
" term2 = np.dot(image_iter.reshape(-1), grad_iter.reshape(-1))\n",
402-
" \n",
393+
"\n",
403394
" term3 = np.dot(grad_iter.reshape(-1), grad_iter.reshape(-1))\n",
404-
" \n",
395+
"\n",
405396
" if niter_lsrtm == 0:\n",
406-
" \n",
397+
"\n",
407398
" alfa = .05 / mmax(grad_full)\n",
408-
" \n",
399+
"\n",
409400
" else:\n",
410401
" abb1 = term1 / term2\n",
411-
" \n",
402+
"\n",
412403
" abb2 = term2 / term3\n",
413-
" \n",
404+
"\n",
414405
" abb3 = abb2 / abb1\n",
415-
" \n",
406+
"\n",
416407
" if abb3 > 0 and abb3 < 1:\n",
417408
" alfa = abb2\n",
418409
" else:\n",
419410
" alfa = abb1\n",
420-
" \n",
421-
" return alfa "
411+
"\n",
412+
" return alfa"
422413
]
423414
},
424415
{
@@ -468,8 +459,8 @@
468459
"\n",
469460
"image = np.zeros((model0.vp.shape[0], model0.vp.shape[1]))\n",
470461
"\n",
471-
"nrec=101\n",
472-
"niter=20 # Number of iterations of the LSRTM\n",
462+
"nrec = 101\n",
463+
"niter = 20 # Number of iterations of the LSRTM\n",
473464
"history = np.zeros((niter, 1)) # Objective function\n",
474465
"\n",
475466
"image_prev = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
@@ -494,7 +485,7 @@
494485
"\n",
495486
" sk = image_up_dev - image_prev\n",
496487
"\n",
497-
" alfa = get_alfa(yk,sk,k)\n",
488+
" alfa = get_alfa(yk, sk, k)\n",
498489
"\n",
499490
" grad_prev = grad_full.data\n",
500491
"\n",
@@ -632,7 +623,7 @@
632623
],
633624
"source": [
634625
"#NBVAL_SKIP\n",
635-
"slices=tuple(slice(model.nbl,-model.nbl) for _ in range(2))\n",
626+
"slices = tuple(slice(model.nbl, -model.nbl) for _ in range(2))\n",
636627
"lsrtm = image_up_dev[slices]\n",
637628
"plot_image(np.diff(lsrtm, axis=1))"
638629
]
@@ -689,13 +680,13 @@
689680
],
690681
"source": [
691682
"#NBVAL_SKIP\n",
692-
"plt.figure(figsize=(8,9))\n",
693-
"x = np.linspace(0,1,101)\n",
694-
"plt.plot(rtm[50,:],x,color=plt.gray(),linewidth=2)\n",
695-
"plt.plot(lsrtm[50,:],x,'r',linewidth=2)\n",
683+
"plt.figure(figsize=(8, 9))\n",
684+
"x = np.linspace(0, 1, 101)\n",
685+
"plt.plot(rtm[50,:], x, color=plt.gray(), linewidth=2)\n",
686+
"plt.plot(lsrtm[50,:],x, 'r',linewidth=2)\n",
696687
"plt.plot(dm_true[50,:],x, 'k--',linewidth=2)\n",
697688
"\n",
698-
"plt.legend(['Initial reflectivity', 'Reflectivity via LSRTM','True Reflectivity'],fontsize=15)\n",
689+
"plt.legend(['Initial reflectivity', 'Reflectivity via LSRTM','True Reflectivity'], fontsize=15)\n",
699690
"plt.ylabel('Depth (Km)')\n",
700691
"plt.xlabel('Amplitude')\n",
701692
"plt.gca().invert_yaxis()\n",
@@ -723,12 +714,12 @@
723714
"source": [
724715
"#NBVAL_SKIP\n",
725716
"time = np.linspace(t0, tn, nt)\n",
726-
"plt.figure(figsize=(8,9))\n",
717+
"plt.figure(figsize=(8, 9))\n",
727718
"plt.ylabel('Time (ms)')\n",
728719
"plt.xlabel('Amplitude')\n",
729-
"plt.plot(d_syn.data[:, 20],time, 'y', label='Calculated data (Last Iteration)')\n",
730-
"plt.plot(d_obs.data[:, 20],time, 'm', label='Observed data')\n",
731-
"plt.legend(loc=\"upper left\",fontsize=12)\n",
720+
"plt.plot(d_syn.data[:, 20], time, 'y', label='Calculated data (Last Iteration)')\n",
721+
"plt.plot(d_obs.data[:, 20], time, 'm', label='Observed data')\n",
722+
"plt.legend(loc=\"upper left\", fontsize=12)\n",
732723
"ax = plt.gca()\n",
733724
"ax.invert_yaxis()\n",
734725
"plt.show()"

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

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -69,9 +69,9 @@
6969
"source": [
7070
"# NBVAL_IGNORE_OUTPUT\n",
7171
"\n",
72-
"shape = (101,101) # 101x101 grid\n",
73-
"spacing = (10.,10.) # spacing of 10 meters\n",
74-
"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",
7575
"nbl = 0 # number of pad points\n",
7676
"\n",
7777
"model = demo_model('layers-tti', spacing=spacing, space_order=8,\n",
@@ -634,7 +634,8 @@
634634
"outputs": [],
635635
"source": [
636636
"# Some useful definitions for plotting if nbl is set to any other value than zero\n",
637-
"nxpad,nzpad = shape[0] + 2 * nbl, shape[1] + 2 * nbl\n",
637+
"nxpad = shape[0] + 2 * nbl,\n",
638+
"nzpad = shape[1] + 2 * nbl\n",
638639
"shape_pad = np.array(shape) + 2 * nbl\n",
639640
"origin_pad = tuple([o - s*nbl for o, s in zip(origin, spacing)])\n",
640641
"extent_pad = tuple([s*(n-1) for s, n in zip(spacing, shape_pad)])"

0 commit comments

Comments
 (0)