Skip to content

Commit

Permalink
examples: Update 13/15/16
Browse files Browse the repository at this point in the history
  • Loading branch information
@georgebisbas
georgebisbas committed Mar 6, 2025
1 parent d92698a commit 19771e2
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 73 deletions.
32 changes: 16 additions & 16 deletions examples/seismic/tutorials/13_LSRTM_acoustic.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -468,41 +468,40 @@
"image = np.zeros((model0.vp.shape[0], model0.vp.shape[1]))\n",
"\n",
"nrec=101\n",
"niter=20 # number of iterations of the LSRTM\n",
"history = np.zeros((niter, 1)) #objective function\n",
"niter=20 # Number of iterations of the LSRTM\n",
"history = np.zeros((niter, 1)) # Objective function\n",
"\n",
"image_prev = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
" \n",
"\n",
"grad_prev = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
"\n",
"yk = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
"\n",
"sk = np.zeros((model0.vp.shape[0],model0.vp.shape[1]))\n",
"\n",
"for k in range(niter) :\n",
" \n",
" dm = image_up_dev # Reflectivity for Calculated data via Born\n",
"for k in range(niter):\n",
"\n",
" dm = image_up_dev # Reflectivity for Calculated data via Born\n",
"\n",
" print('LSRTM Iteration',k+1)\n",
" \n",
"\n",
" objective,grad_full,d_obs,d_syn = lsrtm_gradient(dm)\n",
" \n",
"\n",
" history[k] = objective\n",
" \n",
"\n",
" yk = grad_full.data - grad_prev\n",
" \n",
"\n",
" sk = image_up_dev - image_prev\n",
"\n",
" alfa = get_alfa(yk,sk,k)\n",
" \n",
"\n",
" grad_prev = grad_full.data\n",
"\n",
" image_prev = image_up_dev\n",
" \n",
"\n",
" image_up_dev = image_up_dev - alfa*grad_full.data\n",
" \n",
"\n",
" if k == 0: # Saving the first migration using Born operator.\n",
" \n",
" image = image_up_dev"
]
},
Expand Down Expand Up @@ -559,7 +558,7 @@
" vmin=-.05,\n",
" vmax=.05,\n",
" cmap=cmap,extent=extent)\n",
" \n",
"\n",
" plt.xlabel('X position (km)')\n",
" plt.ylabel('Depth (km)')\n",
"\n",
Expand All @@ -569,6 +568,7 @@
" divider = make_axes_locatable(ax)\n",
" cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
" plt.colorbar(plot, cax=cax)\n",
"\n",
" plt.show()\n"
]
},
Expand Down Expand Up @@ -599,7 +599,7 @@
],
"source": [
"#NBVAL_IGNORE_OUTPUT\n",
"slices=tuple(slice(model.nbl,-model.nbl) for _ in range(2))\n",
"slices=tuple(slice(model.nbl, -model.nbl) for _ in range(2))\n",
"rtm = image[slices]\n",
"plot_image(np.diff(rtm, axis=1))"
]
Expand Down
67 changes: 30 additions & 37 deletions examples/seismic/tutorials/15_tti_qp_pure.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@
"source": [
"import numpy as np\n",
"from devito import (Function, TimeFunction, cos, sin, solve,\n",
" Eq, Operator, configuration, norm)\n",
" Eq, Operator)\n",
"from examples.seismic import TimeAxis, RickerSource, Receiver, demo_model\n",
"from matplotlib import pyplot as plt"
]
Expand All @@ -65,34 +65,23 @@
"execution_count": 2,
"id": "4f545ff1",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Operator `pad_vp` ran in 0.01 s\n",
"Operator `pad_epsilon` ran in 0.01 s\n",
"Operator `pad_delta` ran in 0.01 s\n",
"Operator `pad_theta` ran in 0.01 s\n"
]
}
],
"outputs": [],
"source": [
"# NBVAL_IGNORE_OUTPUT \n",
"# NBVAL_IGNORE_OUTPUT\n",
"\n",
"shape = (101,101) # 101x101 grid\n",
"spacing = (10.,10.) # spacing of 10 meters\n",
"origin = (0.,0.) \n",
"shape = (101,101) # 101x101 grid\n",
"spacing = (10.,10.) # spacing of 10 meters\n",
"origin = (0.,0.)\n",
"nbl = 0 # number of pad points\n",
"\n",
"model = demo_model('layers-tti', spacing=spacing, space_order=8,\n",
" shape=shape, nbl=nbl, nlayers=1)\n",
"\n",
"# initialize Thomsem parameters to those used in Mu et al., (2020)\n",
"model.update('vp', np.ones(shape)*3.6) # km/s\n",
"model.update('vp', np.ones(shape)*3.6) # km/s\n",
"model.update('epsilon', np.ones(shape)*0.23)\n",
"model.update('delta', np.ones(shape)*0.17)\n",
"model.update('theta', np.ones(shape)*(45.*(np.pi/180.))) # radians"
"model.update('theta', np.ones(shape)*(45.*(np.pi/180.))) # radians"
]
},
{
Expand Down Expand Up @@ -233,7 +222,7 @@
"bc += [Eq(pp[t+1,0, z], 0.)]\n",
"bc += [Eq(pp[t+1,shape[0]-1+2*nbl, z], 0.)]\n",
"\n",
"# set source and receivers\n",
"# Set source and receivers\n",
"src = RickerSource(name='src',grid=model.grid,f0=0.02,npoint=1,time_range=time_range)\n",
"src.coordinates.data[:,0] = model.domain_size[0]* .5\n",
"src.coordinates.data[:,1] = model.domain_size[0]* .5\n",
Expand All @@ -243,14 +232,18 @@
"rec = Receiver(name='rec',grid=model.grid,npoint=shape[0],time_range=time_range)\n",
"rec.coordinates.data[:, 0] = np.linspace(model.origin[0],model.domain_size[0], num=model.shape[0])\n",
"rec.coordinates.data[:, 1] = 2*spacing[1]\n",
"\n",
"# Create interpolation expression for receivers\n",
"rec_term = rec.interpolate(expr=p.forward)\n",
"\n",
"# Operators\n",
"optime=Operator([update_p] + src_term + rec_term)\n",
"oppres=Operator([update_q] + bc)\n",
"optime = Operator([update_p] + src_term + rec_term)\n",
"oppres = Operator([update_q] + bc)\n",
"\n",
"\n",
"# you can print the generated code for both operators by typing print(optime) and print(oppres)"
"# You can print the generated code for both operators by uncommenting the following lines\n",
"# print(optime)\n",
"# print(oppres)"
]
},
{
Expand Down Expand Up @@ -620,17 +613,17 @@
],
"source": [
"# NBVAL_IGNORE_OUTPUT\n",
"psave =np.empty ((time_range.num,model.grid.shape[0],model.grid.shape[1]))\n",
"psave =np.empty((time_range.num,model.grid.shape[0], model.grid.shape[1]))\n",
"niter_poisson = 1200\n",
"\n",
"# This is the time loop.\n",
"for step in range(0,time_range.num-2):\n",
" q.data[:,:]=pp.data[(niter_poisson+1)%2,:,:]\n",
"for step in range(0, time_range.num-2):\n",
" q.data[:, :]=pp.data[(niter_poisson+1)%2, :, :]\n",
" optime(time_m=step, time_M=step, dt=dt)\n",
" pp.data[:,:]=0.\n",
" b.data[:,:]=p.data[(step+1)%3,:,:]\n",
" pp.data[:, :]=0.\n",
" b.data[:, :]=p.data[(step+1)%3, :, :]\n",
" oppres(time_M = niter_poisson)\n",
" psave[step,:,:]=p.data[(step+1)%3,:,:]"
" psave[step, :, :]=p.data[(step+1)%3, :, :]"
]
},
{
Expand Down Expand Up @@ -691,18 +684,18 @@
"fig, axes = plt.subplots(2, 5, figsize=(18, 7), sharex=True)\n",
"fig.suptitle(\"Snapshots\", size=14)\n",
"for count, ax in enumerate(axes.ravel()):\n",
" snapshot = factor*count\n",
" snapshot = factor * count\n",
" ax.imshow(np.transpose(psave[snapshot,:,:]), cmap=\"seismic\",\n",
" vmin=-amax, vmax=+amax, extent=plt_extent)\n",
" ax.plot(model.domain_size[0]* .5, model.domain_size[1]* .5, \\\n",
" 'red', linestyle='None', marker='*', markersize=8, label=\"Source\")\n",
" vmin=-amax, vmax=+amax, extent=plt_extent)\n",
" ax.plot(model.domain_size[0] * 0.5, model.domain_size[1] * 0.5,\n",
" 'red', linestyle='None', marker='*', markersize=8, label=\"Source\")\n",
" ax.grid()\n",
" ax.tick_params('both', length=2, width=0.5, which='major',labelsize=10)\n",
" ax.set_title(\"Wavefield at t=%.2fms\" % (factor*count*dt),fontsize=10)\n",
" ax.tick_params('both', length=2, width=0.5, which='major', labelsize=10)\n",
" ax.set_title(f\"Wavefield at t={factor * count * dt:.2f}ms\", fontsize=10)\n",
"for ax in axes[1, :]:\n",
" ax.set_xlabel(\"X Coordinate (m)\",fontsize=10)\n",
" ax.set_xlabel(\"X Coordinate (m)\", fontsize=10)\n",
"for ax in axes[:, 0]:\n",
" ax.set_ylabel(\"Z Coordinate (m)\",fontsize=10)"
" ax.set_ylabel(\"Z Coordinate (m)\", fontsize=10)"
]
},
{
Expand Down
43 changes: 23 additions & 20 deletions examples/seismic/tutorials/16_ader_fd.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,11 +49,12 @@
"source": [
"#NBVAL_IGNORE_OUTPUT\n",
"# Necessary imports\n",
"import devito as dv\n",
"import sympy as sp\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from devito import (Function, TimeFunction, Grid, Operator, Eq, solve,\n",
" VectorTimeFunction, grad, div)\n",
"from examples.seismic import TimeAxis, RickerSource"
]
},
Expand All @@ -73,9 +74,9 @@
"outputs": [],
"source": [
"# 1km x 1km grid\n",
"grid = dv.Grid(shape=(201, 201), extent=(1000., 1000.))\n",
"p = dv.TimeFunction(name='p', grid=grid, space_order=16)\n",
"v = dv.VectorTimeFunction(name='v', grid=grid, space_order=16, staggered=(None, None))"
"grid = Grid(shape=(201, 201), extent=(1000., 1000.))\n",
"p = TimeFunction(name='p', grid=grid, space_order=16)\n",
"v = VectorTimeFunction(name='v', grid=grid, space_order=16, staggered=(None, None))"
]
},
{
Expand All @@ -94,8 +95,8 @@
"outputs": [],
"source": [
"# Material parameters\n",
"c = dv.Function(name='c', grid=grid)\n",
"rho = dv.Function(name='rho', grid=grid)\n",
"c = Function(name='c', grid=grid)\n",
"rho = Function(name='rho', grid=grid)\n",
"\n",
"c.data[:] = 1.5\n",
"c.data[:, :150] = 1.25\n",
Expand Down Expand Up @@ -126,8 +127,8 @@
"metadata": {},
"outputs": [],
"source": [
"# dv.grad(dv.div(v)) is not the same as expanding the continuous operator and then discretising\n",
"# This is because dv.grad(dv.div(v)) applies a gradient stencil to a divergence stencil\n",
"# grad(div(v)) is not the same as expanding the continuous operator and then discretising\n",
"# This is because grad(div(v)) applies a gradient stencil to a divergence stencil\n",
"def graddiv(f):\n",
" return sp.Matrix([[f[0].dx2 + f[1].dxdy],\n",
" [f[0].dxdy + f[1].dy2]])\n",
Expand All @@ -147,8 +148,8 @@
" return f.dx4 + 2*f.dx2dy2 + f.dy4\n",
"\n",
"# First time derivatives\n",
"pdt = rho*c2*dv.div(v)\n",
"vdt = b*dv.grad(p)\n",
"pdt = rho*c2*div(v)\n",
"vdt = b*grad(p)\n",
"\n",
"# Second time derivatives\n",
"pdt2 = c2*p.laplace\n",
Expand Down Expand Up @@ -200,8 +201,8 @@
"dt = grid.stepping_dim.spacing\n",
"\n",
"# Update equations (4th-order ADER timestepping)\n",
"eq_p = dv.Eq(p.forward, p + dt*pdt + (dt**2/2)*pdt2 + (dt**3/6)*pdt3 + (dt**4/24)*pdt4)\n",
"eq_v = dv.Eq(v.forward, v + dt*vdt + (dt**2/2)*vdt2 + (dt**3/6)*vdt3 + (dt**4/24)*vdt4)"
"eq_p = Eq(p.forward, p + dt*pdt + (dt**2/2)*pdt2 + (dt**3/6)*pdt3 + (dt**4/24)*pdt4)\n",
"eq_v = Eq(v.forward, v + dt*vdt + (dt**2/2)*vdt2 + (dt**3/6)*vdt3 + (dt**4/24)*vdt4)"
]
},
{
Expand Down Expand Up @@ -286,7 +287,7 @@
"#NBVAL_IGNORE_OUTPUT\n",
"src_term = src.inject(field=p.forward, expr=src)\n",
"\n",
"op = dv.Operator([eq_p, eq_v] + src_term)\n",
"op = Operator([eq_p, eq_v] + src_term)\n",
"op.apply(dt=op_dt)"
]
},
Expand Down Expand Up @@ -355,20 +356,22 @@
}
],
"source": [
"from devito.types import NODE\n",
"\n",
"#NBVAL_IGNORE_OUTPUT\n",
"ps = dv.TimeFunction(name='ps', grid=grid, space_order=16, staggered=dv.NODE)\n",
"vs = dv.VectorTimeFunction(name='vs', grid=grid, space_order=16)\n",
"ps = TimeFunction(name='ps', grid=grid, space_order=16, staggered=NODE)\n",
"vs = VectorTimeFunction(name='vs', grid=grid, space_order=16)\n",
"\n",
"# First time derivatives\n",
"psdt = rho*c2*dv.div(vs.forward)\n",
"vsdt = b*dv.grad(ps)\n",
"psdt = rho*c2*div(vs.forward)\n",
"vsdt = b*grad(ps)\n",
"\n",
"eq_ps = dv.Eq(ps.forward, ps + dt*psdt)\n",
"eq_vs = dv.Eq(vs.forward, vs + dt*vsdt)\n",
"eq_ps = Eq(ps.forward, ps + dt*psdt)\n",
"eq_vs = Eq(vs.forward, vs + dt*vsdt)\n",
"\n",
"src_term_s = src.inject(field=ps.forward, expr=src)\n",
"\n",
"ops = dv.Operator([eq_vs, eq_ps] + src_term_s)\n",
"ops = Operator([eq_vs, eq_ps] + src_term_s)\n",
"ops.apply(dt=op_dt)"
]
},
Expand Down

0 comments on commit 19771e2

Please sign in to comment.