-
Notifications
You must be signed in to change notification settings - Fork 222
/
elastic_example.py
92 lines (75 loc) · 3.92 KB
/
elastic_example.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
import numpy as np
from argparse import ArgumentParser
from devito import configuration
from devito.logger import info
from examples.seismic.elastic import ElasticWaveSolver
from examples.seismic import demo_model, AcquisitionGeometry
def elastic_setup(shape=(50, 50), spacing=(15.0, 15.0), tn=500., space_order=4,
nbl=10, constant=False, **kwargs):
nrec = 2*shape[0]
preset = 'constant-elastic' if constant else 'layers-elastic'
model = demo_model(preset, space_order=space_order, shape=shape, nbl=nbl,
dtype=kwargs.pop('dtype', np.float32), spacing=spacing)
# Source and receiver geometries
src_coordinates = np.empty((1, len(spacing)))
src_coordinates[0, :] = np.array(model.domain_size) * .5
if len(shape) > 1:
src_coordinates[0, -1] = model.origin[-1] + 2 * spacing[-1]
rec_coordinates = np.empty((nrec, len(spacing)))
rec_coordinates[:, 0] = np.linspace(0., model.domain_size[0], num=nrec)
if len(shape) > 1:
rec_coordinates[:, 1] = np.array(model.domain_size)[1] * .5
rec_coordinates[:, -1] = model.origin[-1] + 2 * spacing[-1]
geometry = AcquisitionGeometry(model, rec_coordinates, src_coordinates,
t0=0.0, tn=tn, src_type='Ricker', f0=0.010)
# Create solver object to provide relevant operators
solver = ElasticWaveSolver(model, geometry, space_order=space_order, **kwargs)
return solver
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])
def test_elastic():
_, _, _, [rec1, rec2, v, tau] = run()
norm = lambda x: np.linalg.norm(x.data.reshape(-1))
assert np.isclose(norm(rec1), 29.6567, atol=1e-3, rtol=0)
assert np.isclose(norm(rec2), 1.19761, atol=1e-3, rtol=0)
if __name__ == "__main__":
description = ("Example script for a set of elastic operators.")
parser = ArgumentParser(description=description)
parser.add_argument('--2d', dest='dim2', default=False, action='store_true',
help="Preset to determine the physical problem setup")
parser.add_argument('-a', '--autotune', default='off',
choices=(configuration._accepted['autotuning']),
help="Operator auto-tuning mode")
parser.add_argument("-so", "--space_order", default=4,
type=int, help="Space order of the simulation")
parser.add_argument("--nbl", default=40,
type=int, help="Number of boundary layers around the domain")
parser.add_argument("-dse", default="advanced",
choices=["noop", "basic", "advanced", "aggressive"],
help="Devito symbolic engine (DSE) mode")
parser.add_argument("-dle", default="advanced",
choices=["noop", "advanced", "speculative"],
help="Devito loop engine (DLEE) mode")
parser.add_argument("--constant", default=False, action='store_true',
help="Constant velocity model, default is a two layer model")
args = parser.parse_args()
# 2D preset parameters
if args.dim2:
shape = (150, 150)
spacing = (10.0, 10.0)
tn = 750.0
# 3D preset parameters
else:
shape = (150, 150, 150)
spacing = (10.0, 10.0, 10.0)
tn = 1250.0
run(shape=shape, spacing=spacing, nbl=args.nbl, tn=tn, dle=args.dle,
space_order=args.space_order, autotune=args.autotune, constant=args.constant,
dse=args.dse)