Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# Compute the full wavefield u0
_, u0, _ = solver.forward(save=True, vp=model0.vp)
# Compute initial born perturbation from m - m0
dm = (solver.model.vp.data**(-2) - model0.vp.data**(-2))
du, _, _, _ = solver.born(dm, vp=model0.vp)
# Compute gradientfrom initial perturbation
im, _ = solver.gradient(du, u0, vp=model0.vp)
# Adjoint test: Verify matches closely
term1 = np.dot(im.data.reshape(-1), dm.reshape(-1))
term2 = linalg.norm(du.data.reshape(-1))**2
info(': %f, : %f, difference: %4.4e, ratio: %f'
% (term1, term2, (term1 - term2)/term1, term1 / term2))
assert np.isclose((term1 - term2)/term1, 0., atol=1.e-12)
nbpml=nbpml, tn=tn, space_order=space_order,
**(presets[mkey]), dtype=np.float64)
# Create adjoint receiver symbol
srca = Receiver(name='srca', grid=solver.model.grid,
time_range=solver.geometry.time_axis,
coordinates=solver.geometry.src_positions)
# Run forward and adjoint operators
rec, _, _ = solver.forward(save=False)
solver.adjoint(rec=rec, srca=srca)
# Adjoint test: Verify matches closely
term1 = np.dot(srca.data.reshape(-1), solver.geometry.src.data)
term2 = linalg.norm(rec.data.reshape(-1)) ** 2
info(': %f, : %f, difference: %4.4e, ratio: %f'
% (term1, term2, (term1 - term2)/term1, term1 / term2))
assert np.isclose((term1 - term2)/term1, 0., atol=1.e-12)
def run(shape=(50, 50), spacing=(20.0, 20.0), tn=1000.0,
space_order=4, nbl=40, autotune=False, constant=False, **kwargs):
solver = elastic_setup(shape=shape, spacing=spacing, nbl=nbl, tn=tn,
space_order=space_order, constant=constant, **kwargs)
info("Applying Forward")
# Define receiver geometry (spread across x, just below surface)
rec1, rec2, v, tau, summary = solver.forward(autotune=autotune)
return (summary.gflopss, summary.oi, summary.timings,
[rec1, rec2, v, tau])
if self.profile:
args.append(self.profiler.as_ctypes_pointer(Profiler.TIME))
if isinstance(self.compiler, IntelMICCompiler):
# Off-load propagator kernel via pymic stream
self.compiler._stream.invoke(f, *args)
self.compiler._stream.sync()
else:
f(*args)
if self.profile:
if verbose:
shape = str(self.shape).replace(', ', ' x ')
cb = str(self.block_sizes) if self.cache_blocking else 'None'
info("Shape: %s - Cache Blocking: %s" % (shape, cb))
info("Time: %f s (%s MCells/s)" % (self.total_time, self.mcells))
key = LOOP_BODY.name
info("Stencil: %f OI, %.2f GFlops/s (time: %f s)" %
(self.oi[key], self.gflopss[key], self.timings[key]))
if preset == 'constant':
# With a new m as Constant
v0 = Constant(name="v", value=2.0, dtype=np.float32)
solver.forward(save=save, vp=v0)
# With a new vp as a scalar value
solver.forward(save=save, vp=2.0)
if not full_run:
return summary.gflopss, summary.oi, summary.timings, [rec, u.data]
# Smooth velocity
initial_vp = Function(name='v0', grid=solver.model.grid, space_order=space_order)
smooth(initial_vp, solver.model.vp)
dm = np.float32(initial_vp.data**(-2) - solver.model.vp.data**(-2))
info("Applying Adjoint")
solver.adjoint(rec, autotune=autotune)
info("Applying Born")
solver.born(dm, autotune=autotune)
info("Applying Gradient")
solver.gradient(rec, u, autotune=autotune, checkpointing=checkpointing)
return summary.gflopss, summary.oi, summary.timings, [rec, u.data]
source_mask = Function(name='source_mask', shape=shape, dimensions=(x, y, z),
dtype=np.int32)
source_id = Function(name='source_id', grid=model.grid, dtype=np.int32)
source_id.data[nzinds[0], nzinds[1], nzinds[2]] = tuple(np.arange(1, len(nzinds[0])+1))
source_mask.data[nzinds[0], nzinds[1], nzinds[2]] = 1
plot3d(source_mask.data, model)
print("Number of unique affected points is:", len(nzinds[0]))
assert(source_id.data[nzinds[0][0], nzinds[1][0], nzinds[2][0]] == 1)
assert(source_id.data[nzinds[0][-1], nzinds[1][-1], nzinds[2][-1]] == len(nzinds[0]))
assert(source_id.data[nzinds[0][len(nzinds[0])-1], nzinds[1][len(nzinds[0])-1],
nzinds[2][len(nzinds[0])-1]] == len(nzinds[0]))
info("---Source_mask and source_id is built here-------")
nnz_shape = (model.grid.shape[0], model.grid.shape[1]) # Change only 3rd dim
nnz_sp_source_mask = TimeFunction(name='nnz_sp_source_mask',
shape=(list(shape[:2])),
dimensions=(x, y), dtype=np.int32)
nnz_sp_source_mask.data[:, :] = source_mask.data.sum(2)
inds = np.where(source_mask.data == 1)
maxz = len(np.unique(inds[2]))
sparse_shape = (model.grid.shape[0], model.grid.shape[1], maxz) # Change only 3rd dim
assert(len(nnz_sp_source_mask.dimensions) == 2)
# Note:sparse_source_id is not needed as long as sparse info is kept in mask
# sp_source_id.data[inds[0],inds[1],:] = inds[2][:maxz]
source_mask = Function(name='source_mask', shape=shape, dimensions=(x, y, z),
dtype=np.int32)
source_id = Function(name='source_id', grid=model.grid, dtype=np.int32)
source_id.data[nzinds[0], nzinds[1], nzinds[2]] = tuple(np.arange(1, len(nzinds[0])+1))
source_mask.data[nzinds[0], nzinds[1], nzinds[2]] = 1
plot3d(source_mask.data, model)
print("Number of unique affected points is:", len(nzinds[0]))
assert(source_id.data[nzinds[0][0], nzinds[1][0], nzinds[2][0]] == 1)
assert(source_id.data[nzinds[0][-1], nzinds[1][-1], nzinds[2][-1]] == len(nzinds[0]))
assert(source_id.data[nzinds[0][len(nzinds[0])-1], nzinds[1][len(nzinds[0])-1],
nzinds[2][len(nzinds[0])-1]] == len(nzinds[0]))
info("---Source_mask and source_id is built here-------")
nnz_shape = (model.grid.shape[0], model.grid.shape[1]) # Change only 3rd dim
nnz_sp_source_mask = TimeFunction(name='nnz_sp_source_mask',
shape=([1] + list(shape[:2])),
dimensions=(time, x, y), time_order=0, dtype=np.int32)
nnz_sp_source_mask.data[0, :, :] = source_mask.data.sum(2)
inds = np.where(source_mask.data == 1)
maxz = len(np.unique(inds[2]))
sparse_shape = (model.grid.shape[0], model.grid.shape[1], maxz) # Change only 3rd dim
assert(len(nnz_sp_source_mask.dimensions) == 3)
sp_source_mask = TimeFunction(name='sp_source_mask', shape=([1] + list(sparse_shape)),
dimensions=(time, x, y, z), time_order=0, dtype=np.int32)
def run(shape=(50, 50), spacing=(20.0, 20.0), tn=1000.0,
space_order=4, nbl=40, autotune=False, constant=False, **kwargs):
solver = viscoelastic_setup(shape=shape, spacing=spacing, nbl=nbl, tn=tn,
space_order=space_order, constant=constant, **kwargs)
info("Applying Forward")
# Define receiver geometry (spread across x, just below surface)
rec1, rec2, v, tau, summary = solver.forward(autotune=autotune)
return (summary.gflopss, summary.oi, summary.timings,
[rec1, rec2, v, tau])
def run(shape=(50, 50), spacing=(20.0, 20.0), tn=1000.0,
space_order=4, nbl=40, autotune=False, constant=False, **kwargs):
solver = elastic_setup(shape=shape, spacing=spacing, nbl=nbl, tn=tn,
space_order=space_order, constant=constant, **kwargs)
info("Applying Forward")
# Define receiver geometry (spread across x, just below surface)
rec1, rec2, v, tau, summary = solver.forward(autotune=autotune)
return (summary.gflopss, summary.oi, summary.timings,
[rec1, rec2, v, tau])
source_mask = Function(name='source_mask', shape=shape, dimensions=(x, y, z),
dtype=np.int32)
source_id = Function(name='source_id', grid=model.grid, dtype=np.int32)
source_id.data[nzinds[0], nzinds[1], nzinds[2]] = tuple(np.arange(1, len(nzinds[0])+1))
source_mask.data[nzinds[0], nzinds[1], nzinds[2]] = 1
# plot3d(source_mask.data, model)
print("Number of unique affected points is:", len(nzinds[0]))
assert(source_id.data[nzinds[0][0], nzinds[1][0], nzinds[2][0]] == 1)
assert(source_id.data[nzinds[0][-1], nzinds[1][-1], nzinds[2][-1]] == len(nzinds[0]))
assert(source_id.data[nzinds[0][len(nzinds[0])-1], nzinds[1][len(nzinds[0])-1],
nzinds[2][len(nzinds[0])-1]] == len(nzinds[0]))
info("---Source_mask and source_id is built here-------")
nnz_shape = (model.grid.shape[0], model.grid.shape[1]) # Change only 3rd dim
nnz_sp_source_mask = TimeFunction(name='nnz_sp_source_mask',
shape=(list(shape[:2])),
dimensions=(x, y), dtype=np.int32)
nnz_sp_source_mask.data[:, :] = source_mask.data.sum(2)
inds = np.where(source_mask.data == 1)
maxz = len(np.unique(inds[2]))
sparse_shape = (model.grid.shape[0], model.grid.shape[1], maxz) # Change only 3rd dim
assert(len(nnz_sp_source_mask.dimensions) == 2)
# Note:sparse_source_id is not needed as long as sparse info is kept in mask
# sp_source_id.data[inds[0],inds[1],:] = inds[2][:maxz]