Skip to content

Commit 1091ee7

Browse files
author
George Bisbas
committed
examples: Further reviews and updates
1 parent a844834 commit 1091ee7

2 files changed

Lines changed: 51 additions & 59 deletions

File tree

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

Lines changed: 46 additions & 55 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,56 +300,56 @@
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
" for i in range(nshots):\n",
333-
" \n",
324+
"\n",
334325
" #Observed Data using Born's operator\n",
335326
" geometry.src_positions[0, :] = source_locations[i, :]\n",
336327
"\n",
337328
" _, u0, _ = solver.forward(vp=model0.vp, save=True)\n",
338-
" \n",
329+
"\n",
339330
" _, _, _,_ = solver.jacobian(dm_true, vp=model0.vp, rec = d_obs)\n",
340-
" \n",
331+
"\n",
341332
" #Calculated Data using Born's operator\n",
342333
" solver.jacobian(dm, vp=model0.vp, rec = d_syn)\n",
343-
" \n",
334+
"\n",
344335
" residual.data[:] = d_syn.data[:]- d_obs.data[:]\n",
345-
" \n",
336+
"\n",
346337
" grad_shot,_ = solver.gradient(rec=residual, u=u0, vp=model0.vp)\n",
347-
" \n",
338+
"\n",
348339
" src_illum_upd = Eq(src_illum, src_illum + u0**2)\n",
349340
" op_src = Operator([src_illum_upd])\n",
350341
" op_src.apply()\n",
351-
" \n",
342+
"\n",
352343
" grad_sum = Eq(grad_full, grad_full + grad_shot)\n",
353344
" op_grad = Operator([grad_sum])\n",
354345
" op_grad.apply()\n",
355-
" \n",
346+
"\n",
356347
" objective += .5*norm(residual)**2\n",
357-
" \n",
348+
"\n",
358349
" grad_f = Eq(grad_illum, grad_full/(src_illum+10**-9))\n",
359350
" op_gradf = Operator([grad_f])\n",
360351
" op_gradf.apply()\n",
361-
" \n",
352+
"\n",
362353
" return objective,grad_illum,d_obs,d_syn"
363354
]
364355
},
@@ -393,31 +384,31 @@
393384
"outputs": [],
394385
"source": [
395386
"def get_alfa(grad_iter,image_iter,niter_lsrtm):\n",
396-
" \n",
397-
" \n",
387+
"\n",
388+
"\n",
398389
" term1 = np.dot(image_iter.reshape(-1), image_iter.reshape(-1))\n",
399-
" \n",
390+
"\n",
400391
" term2 = np.dot(image_iter.reshape(-1), grad_iter.reshape(-1))\n",
401-
" \n",
392+
"\n",
402393
" term3 = np.dot(grad_iter.reshape(-1), grad_iter.reshape(-1))\n",
403-
" \n",
394+
"\n",
404395
" if niter_lsrtm == 0:\n",
405-
" \n",
396+
"\n",
406397
" alfa = .05 / mmax(grad_full)\n",
407-
" \n",
398+
"\n",
408399
" else:\n",
409400
" abb1 = term1 / term2\n",
410-
" \n",
401+
"\n",
411402
" abb2 = term2 / term3\n",
412-
" \n",
403+
"\n",
413404
" abb3 = abb2 / abb1\n",
414-
" \n",
405+
"\n",
415406
" if abb3 > 0 and abb3 < 1:\n",
416407
" alfa = abb2\n",
417408
" else:\n",
418409
" alfa = abb1\n",
419-
" \n",
420-
" return alfa "
410+
"\n",
411+
" return alfa"
421412
]
422413
},
423414
{
@@ -467,8 +458,8 @@
467458
"\n",
468459
"image = np.zeros((model0.vp.shape[0], model0.vp.shape[1]))\n",
469460
"\n",
470-
"nrec=101\n",
471-
"niter=20 # Number of iterations of the LSRTM\n",
461+
"nrec = 101\n",
462+
"niter = 20 # Number of iterations of the LSRTM\n",
472463
"history = np.zeros((niter, 1)) # Objective function\n",
473464
"\n",
474465
"image_prev = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
@@ -493,7 +484,7 @@
493484
"\n",
494485
" sk = image_up_dev - image_prev\n",
495486
"\n",
496-
" alfa = get_alfa(yk,sk,k)\n",
487+
" alfa = get_alfa(yk, sk, k)\n",
497488
"\n",
498489
" grad_prev = grad_full.data\n",
499490
"\n",
@@ -631,7 +622,7 @@
631622
],
632623
"source": [
633624
"#NBVAL_SKIP\n",
634-
"slices=tuple(slice(model.nbl,-model.nbl) for _ in range(2))\n",
625+
"slices = tuple(slice(model.nbl, -model.nbl) for _ in range(2))\n",
635626
"lsrtm = image_up_dev[slices]\n",
636627
"plot_image(np.diff(lsrtm, axis=1))"
637628
]
@@ -688,13 +679,13 @@
688679
],
689680
"source": [
690681
"#NBVAL_SKIP\n",
691-
"plt.figure(figsize=(8,9))\n",
692-
"x = np.linspace(0,1,101)\n",
693-
"plt.plot(rtm[50,:],x,color=plt.gray(),linewidth=2)\n",
694-
"plt.plot(lsrtm[50,:],x,'r',linewidth=2)\n",
682+
"plt.figure(figsize=(8, 9))\n",
683+
"x = np.linspace(0, 1, 101)\n",
684+
"plt.plot(rtm[50,:], x, color=plt.gray(), linewidth=2)\n",
685+
"plt.plot(lsrtm[50,:],x, 'r',linewidth=2)\n",
695686
"plt.plot(dm_true[50,:],x, 'k--',linewidth=2)\n",
696687
"\n",
697-
"plt.legend(['Initial reflectivity', 'Reflectivity via LSRTM','True Reflectivity'],fontsize=15)\n",
688+
"plt.legend(['Initial reflectivity', 'Reflectivity via LSRTM','True Reflectivity'], fontsize=15)\n",
698689
"plt.ylabel('Depth (Km)')\n",
699690
"plt.xlabel('Amplitude')\n",
700691
"plt.gca().invert_yaxis()\n",
@@ -722,12 +713,12 @@
722713
"source": [
723714
"#NBVAL_SKIP\n",
724715
"time = np.linspace(t0, tn, nt)\n",
725-
"plt.figure(figsize=(8,9))\n",
716+
"plt.figure(figsize=(8, 9))\n",
726717
"plt.ylabel('Time (ms)')\n",
727718
"plt.xlabel('Amplitude')\n",
728-
"plt.plot(d_syn.data[:, 20],time, 'y', label='Calculated data (Last Iteration)')\n",
729-
"plt.plot(d_obs.data[:, 20],time, 'm', label='Observed data')\n",
730-
"plt.legend(loc=\"upper left\",fontsize=12)\n",
719+
"plt.plot(d_syn.data[:, 20], time, 'y', label='Calculated data (Last Iteration)')\n",
720+
"plt.plot(d_obs.data[:, 20], time, 'm', label='Observed data')\n",
721+
"plt.legend(loc=\"upper left\", fontsize=12)\n",
731722
"ax = plt.gca()\n",
732723
"ax.invert_yaxis()\n",
733724
"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)