/
wavesolver.py
128 lines (112 loc) · 5.16 KB
/
wavesolver.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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
# coding: utf-8
from devito import TimeFunction, warning
from devito.tools import memoized_meth
from examples.seismic.tti.operators import ForwardOperator, particle_velocity_fields
from examples.seismic import Receiver
class AnisotropicWaveSolver(object):
"""
Solver object that provides operators for seismic inversion problems
and encapsulates the time and space discretization for a given problem
setup.
Parameters
----------
model : Model
Object containing the physical parameters.
geometry : AcquisitionGeometry
Geometry object that contains the source (SparseTimeFunction) and
receivers (SparseTimeFunction) and their position.
space_order : int, optional
Order of the spatial stencil discretisation. Defaults to 4.
Notes
-----
space_order must be even and it is recommended to be a multiple of 4
"""
def __init__(self, model, geometry, space_order=4, **kwargs):
self.model = model
self.geometry = geometry
if space_order % 2 != 0:
raise ValueError("space_order must be even but got %s" % space_order)
if space_order % 4 != 0:
warning("It is recommended for space_order to be a multiple of 4 " +
"but got %s" % space_order)
self.space_order = space_order
self.dt = self.model.critical_dt
# Cache compiler options
self._kwargs = kwargs
@memoized_meth
def op_fwd(self, kernel='centered', save=False):
"""Cached operator for forward runs with buffered wavefield"""
return ForwardOperator(self.model, save=save, geometry=self.geometry,
space_order=self.space_order,
kernel=kernel, **self._kwargs)
def forward(self, src=None, rec=None, u=None, v=None, vp=None,
epsilon=None, delta=None, theta=None, phi=None,
save=False, kernel='centered', **kwargs):
"""
Forward modelling function that creates the necessary
data objects for running a forward modelling operator.
Parameters
----------
geometry : AcquisitionGeometry
Geometry object that contains the source (SparseTimeFunction) and
receivers (SparseTimeFunction) and their position.
u : TimeFunction, optional
The computed wavefield first component.
v : TimeFunction, optional
The computed wavefield second component.
vp : Function or float, optional
The time-constant velocity.
epsilon : Function or float, optional
The time-constant first Thomsen parameter.
delta : Function or float, optional
The time-constant second Thomsen parameter.
theta : Function or float, optional
The time-constant Dip angle (radians).
phi : Function or float, optional
The time-constant Azimuth angle (radians).
save : int or Buffer
Option to store the entire (unrolled) wavefield.
kernel : str, optional
Type of discretization, centered or shifted.
Returns
-------
Receiver, wavefield and performance summary.
"""
if kernel == 'staggered':
time_order = 1
dims = self.model.space_dimensions
stagg_u = (-dims[-1])
stagg_v = (-dims[0], -dims[1]) if self.model.grid.dim == 3 else (-dims[0])
else:
time_order = 2
stagg_u = stagg_v = None
# Source term is read-only, so re-use the default
src = src or self.geometry.src
# Create a new receiver object to store the result
rec = rec or Receiver(name='rec', grid=self.model.grid,
time_range=self.geometry.time_axis,
coordinates=self.geometry.rec_positions)
# Create the forward wavefield if not provided
if u is None:
u = TimeFunction(name='u', grid=self.model.grid, staggered=stagg_u,
save=self.geometry.nt if save else None,
time_order=time_order, space_order=self.space_order)
# Create the forward wavefield if not provided
if v is None:
v = TimeFunction(name='v', grid=self.model.grid, staggered=stagg_v,
save=self.geometry.nt if save else None,
time_order=time_order, space_order=self.space_order)
if kernel == 'staggered':
vx, vz, vy = particle_velocity_fields(self.model, self.space_order)
kwargs["vx"] = vx
kwargs["vz"] = vz
if vy is not None:
kwargs["vy"] = vy
# Pick vp and Thomsen parameters from model unless explicitly provided
kwargs.update(self.model.physical_params(vp=vp, epsilon=epsilon, delta=delta,
theta=theta, phi=phi))
# Execute operator and return wavefield and receiver data
op = self.op_fwd(kernel, save)
summary = op.apply(src=src, rec=rec, u=u, v=v,
dt=kwargs.pop('dt', self.dt), **kwargs)
return rec, u, v, summary