Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gf/tractions.py: 70%
166 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-06 06:59 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Traction field parameterizations used by :doc:`PDR source
8</topics/pseudo-dynamic-rupture>`.
9'''
11import logging
12import numpy as num
13from pyrocko.guts import Object, Float, List, StringChoice, Int
14from pyrocko.guts_array import Array
16logger = logging.getLogger('pyrocko.gf.tractions')
17km = 1e3
18d2r = num.pi/180.
19r2d = 180./num.pi
22def tukey_window(N, alpha):
23 assert alpha <= 1.
24 window = num.ones(N)
25 n = num.arange(N)
27 N_f = int((alpha * N)//2)
28 window[:N_f] = .5 * (1. - num.cos((2*num.pi * n[:N_f])/(alpha * N)))
29 window[(N-N_f):] = window[:N_f][::-1]
30 return window
33def planck_window(N, epsilon):
34 assert epsilon <= 1.
35 window = num.ones(N)
36 n = num.arange(N)
38 N_f = int((epsilon * N))
39 window[:N_f] = \
40 (1. + num.exp((epsilon * N) / n[:N_f] -
41 (epsilon * N) / ((epsilon * N - n[:N_f]))))**-1.
42 window[(N-N_f):] = window[:N_f][::-1]
43 return window
46class AbstractTractionField(Object):
47 '''
48 Base class for multiplicative traction fields (tapers).
50 Fields of this type a re multiplied in the
51 :py:class:`~pyrocko.gf.tractions.TractionComposition`
52 '''
53 operation = 'mult'
55 def get_tractions(self, nx, ny, patches):
56 raise NotImplementedError
59class TractionField(AbstractTractionField):
60 '''
61 Base class for additive traction fields.
63 Fields of this type are added in the
64 :py:class:`~pyrocko.gf.tractions.TractionComposition`
65 '''
66 operation = 'add'
68 def get_tractions(self, nx, ny, patches):
69 raise NotImplementedError
72class TractionComposition(TractionField):
73 '''
74 Composition of traction fields.
76 :py:class:`~pyrocko.gf.tractions.TractionField` and
77 :py:class:`~pyrocko.gf.tractions.AbstractTractionField` can be combined
78 to realize a combination of different fields.
79 '''
80 components = List.T(
81 AbstractTractionField.T(),
82 default=[],
83 help='Ordered list of tractions.')
85 def get_tractions(self, nx, ny, patches=None):
86 npatches = nx * ny
87 tractions = num.zeros((npatches, 3))
89 for comp in self.components:
90 if comp.operation == 'add':
91 tractions += comp.get_tractions(nx, ny, patches)
92 elif comp.operation == 'mult':
93 tractions *= comp.get_tractions(nx, ny, patches)
94 else:
95 raise AttributeError(
96 'Component %s has an invalid operation %s.' %
97 (comp, comp.operation))
99 return tractions
101 def add_component(self, field):
102 logger.debug('Adding traction component.')
103 self.components.append(field)
106class HomogeneousTractions(TractionField):
107 '''
108 Homogeneous traction field.
110 The traction vectors in strike, dip and normal direction are acting
111 homogeneously on the rupture plane.
112 '''
114 strike = Float.T(
115 default=1.,
116 help='Tractions in strike direction [Pa].')
117 dip = Float.T(
118 default=1.,
119 help='Traction in dip direction (up) [Pa].')
120 normal = Float.T(
121 default=1.,
122 help='Traction in normal direction [Pa].')
124 def get_tractions(self, nx, ny, patches=None):
125 npatches = nx * ny
127 return num.tile(
128 (self.strike, self.dip, self.normal), npatches) \
129 .reshape(-1, 3)
132class DirectedTractions(TractionField):
133 '''
134 Directed traction field.
136 The traction vectors are following a uniform ``rake``.
137 '''
139 rake = Float.T(
140 default=0.,
141 help='Rake angle in [deg], '
142 'measured counter-clockwise from right-horizontal '
143 'in on-plane view. Rake is translated into homogenous tractions '
144 'in strike and up-dip direction.')
145 traction = Float.T(
146 default=1.,
147 help='Traction in rake direction [Pa].')
149 def get_tractions(self, nx, ny, patches=None):
150 npatches = nx * ny
152 strike = num.cos(self.rake*d2r) * self.traction
153 dip = num.sin(self.rake*d2r) * self.traction
154 normal = 0.
156 return num.tile((strike, dip, normal), npatches).reshape(-1, 3)
159class SelfSimilarTractions(TractionField):
160 '''
161 Traction model following Power & Tullis (1991).
163 The traction vectors are calculated as a sum of 2D-cosines with a constant
164 amplitude / wavelength ratio. The wavenumber kx and ky are constant for
165 each cosine function. The rank defines the maximum wavenumber used for
166 summation. So, e.g. a rank of 3 will lead to a summation of cosines with
167 ``kx = ky`` in (1, 2, 3).
168 Each cosine has an associated phases, which defines both the phase shift
169 and also the shift from the rupture plane centre.
170 Finally the summed cosines are translated into shear tractions based on the
171 rake and normalized with ``traction_max``.
173 '''
174 rank = Int.T(
175 default=1,
176 help='Maximum summed cosine wavenumber/spatial frequency.')
178 rake = Float.T(
179 default=0.,
180 help='Rake angle in [deg], '
181 'measured counter-clockwise from right-horizontal '
182 'in on-plane view. Rake is translated into homogenous tractions '
183 'in strike and up-dip direction.')
185 traction_max = Float.T(
186 default=1.,
187 help='Maximum traction vector length [Pa].')
189 phases = Array.T(
190 optional=True,
191 dtype=num.float64,
192 shape=(None,),
193 help='Phase shift of the cosines in [rad].')
195 def get_phases(self):
196 if self.phases is not None:
197 if self.phases.shape[0] == self.rank:
198 return self.phases
200 return (num.random.random(self.rank) * 2. - 1.) * num.pi
202 def get_tractions(self, nx, ny, patches=None):
203 z = num.zeros((ny, nx))
204 phases = self.get_phases()
206 for i in range(1, self.rank+1):
207 x = num.linspace(-i*num.pi, i*num.pi, nx) + i*phases[i-1]
208 y = num.linspace(-i*num.pi, i*num.pi, ny) + i*phases[i-1]
209 x, y = num.meshgrid(x, y)
210 r = num.sqrt(x**2 + y**2)
211 z += 1. / i * num.cos(r + phases[i-1])
213 t = num.zeros((nx*ny, 3))
214 t[:, 0] = num.cos(self.rake*d2r) * z.ravel(order='F')
215 t[:, 1] = num.sin(self.rake*d2r) * z.ravel(order='F')
217 t *= self.traction_max / num.max(num.linalg.norm(t, axis=1))
219 return t
222class FractalTractions(TractionField):
223 '''
224 Fractal traction field.
225 '''
227 rseed = Int.T(
228 default=None,
229 optional=True,
230 help='Seed for :py:class:`~numpy.random.RandomState`.'
231 'If ``None``, an random seed will be initialized.')
233 rake = Float.T(
234 default=0.,
235 help='Rake angle in [deg], '
236 'measured counter-clockwise from right-horizontal '
237 'in on-plane view. Rake is translated into homogenous tractions '
238 'in strike and up-dip direction.')
240 traction_max = Float.T(
241 default=1.,
242 help='Maximum traction vector length [Pa].')
244 def __init__(self, *args, **kwargs):
245 super().__init__(*args, **kwargs)
246 if self.rseed is None:
247 self.rseed = num.random.randint(0, 2**31-1)
248 self._data = None
250 def _get_data(self, nx, ny):
251 if self._data is None:
252 rstate = num.random.RandomState(self.rseed)
253 self._data = rstate.rand(nx, ny)
255 return self._data
257 def get_tractions(self, nx, ny, patches=None):
258 if patches is None:
259 raise AttributeError(
260 'Patches needs to be given for this traction field.')
261 npatches = nx * ny
262 dx = -patches[0].al1 + patches[0].al2
263 dy = -patches[0].aw1 + patches[0].aw2
265 # Create random data and get spectrum and power spectrum
266 data = self._get_data(nx, ny)
267 spec = num.fft.fftshift(num.fft.fft2(data))
268 power_spec = (num.abs(spec)/spec.size)**2
270 # Get 0-centered wavenumbers (k_rad == 0.) is in the centre
271 kx = num.fft.fftshift(num.fft.fftfreq(nx, d=dx))
272 ky = num.fft.fftshift(num.fft.fftfreq(ny, d=dy))
273 k_rad = num.sqrt(ky[:, num.newaxis]**2 + kx[num.newaxis, :]**2)
275 # Define wavenumber bins
276 k_bins = num.arange(0, num.max(k_rad), num.max(k_rad)/10.)
278 # Set amplitudes within wavenumber bins to power_spec * 1 / k_max
279 amps = num.zeros(k_rad.shape)
280 amps[k_rad == 0.] = 1.
282 for i in range(k_bins.size-1):
283 k_min = k_bins[i]
284 k_max = k_bins[i+1]
285 r = num.logical_and(k_rad > k_min, k_rad <= k_max)
286 amps[r] = power_spec.T[r]
287 amps = num.sqrt(amps * data.size * num.pi * 4)
289 amps[k_rad > k_bins.max()] = power_spec.ravel()[num.argmax(power_spec)]
291 # Multiply spectrum by amplitudes and inverse fft into demeaned noise
292 spec *= amps.T
294 tractions = num.abs(num.fft.ifft2(spec))
295 tractions -= num.mean(tractions)
296 tractions *= self.traction_max / num.abs(tractions).max()
298 t = num.zeros((npatches, 3))
299 t[:, 0] = num.cos(self.rake*d2r) * tractions.ravel(order='C')
300 t[:, 1] = num.sin(self.rake*d2r) * tractions.ravel(order='C')
302 return t
305class RectangularTaper(AbstractTractionField):
306 width = Float.T(
307 default=.2,
308 help='Width of the taper as a fraction of the plane.')
310 type = StringChoice.T(
311 choices=('tukey', ),
312 default='tukey',
313 help='Type of the taper, default: "tukey".')
315 def get_tractions(self, nx, ny, patches=None):
316 if self.type == 'tukey':
317 x = tukey_window(nx, self.width)
318 y = tukey_window(ny, self.width)
319 return (x[:, num.newaxis] * y).ravel()[:, num.newaxis]
321 raise AttributeError('Unknown type: %s' % self.type)
324class DepthTaper(AbstractTractionField):
325 depth_start = Float.T(
326 help='Depth where the taper begins [m].')
328 depth_stop = Float.T(
329 help='Depth where taper ends and drops to zero [m].')
331 type = StringChoice.T(
332 choices=('linear', ),
333 default='linear',
334 help='Type of the taper, default: "linear".')
336 def get_tractions(self, nx, ny, patches):
337 assert self.depth_stop > self.depth_start
338 depths = num.array([p.depth for p in patches])
340 if self.type == 'linear':
341 slope = self.depth_stop - self.depth_start
342 depths -= self.depth_stop
343 depths /= -slope
344 depths[depths > 1.] = 1.
345 depths[depths < 0.] = 0.
346 return depths[:, num.newaxis]
349def plot_tractions(tractions, nx=15, ny=12, depth=10*km, component='strike'):
350 '''
351 Plot traction model for quick inspection.
353 :param tractions:
354 Traction field or traction composition to be displayed.
355 :type tractions:
356 :py:class:`pyrocko.gf.tractions.TractionField`
358 :param nx:
359 Number of patches along strike.
360 :type nx:
361 int
363 :param ny:
364 Number of patches down dip.
365 :type ny:
366 int
368 :param depth:
369 Depth of the rupture plane center in [m].
370 :type depth:
371 float
373 :param component:
374 Choice of traction component to be shown. Available: ``'tx'`` (along
375 strike), ``'ty'`` (up dip), ``'tz'`` (normal), ``'absolute'`` (vector
376 length).
377 :type component:
378 str
379 '''
380 import matplotlib.pyplot as plt
381 from pyrocko.modelling.okada import OkadaSource
383 comp2idx = dict(
384 tx=0, ty=1, tz=2)
386 source = OkadaSource(
387 lat=0.,
388 lon=0.,
389 depth=depth,
390 al1=-20*km, al2=20*km,
391 aw1=-15*km, aw2=15*km,
392 strike=120., dip=90., rake=90.,
393 slip=5.)
395 patches, _ = source.discretize(nx, ny)
396 tractions = tractions.get_tractions(nx, ny, patches)
398 if component in comp2idx:
399 tractions = tractions[:, comp2idx[component]].reshape(nx, ny)
400 elif component == 'absolute':
401 tractions = num.linalg.norm(tractions, axis=1).reshape(nx, ny)
402 else:
403 raise ValueError('Given component is not valid.')
405 fig = plt.figure()
406 ax = fig.gca()
408 ax.imshow(tractions)
410 plt.show()
413__all__ = [
414 'AbstractTractionField',
415 'TractionField',
416 'TractionComposition',
417 'HomogeneousTractions',
418 'DirectedTractions',
419 'FractalTractions',
420 'SelfSimilarTractions',
421 'RectangularTaper',
422 'DepthTaper',
423 'plot_tractions']