1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
7import logging
9from pyrocko import moment_tensor as mt
10from pyrocko.guts import Float, String, Timestamp, Int
11from pyrocko.model import Location
12from pyrocko.modelling import okada_ext
13from pyrocko.util import get_threadpool_limits
15guts_prefix = 'modelling'
17logger = logging.getLogger(__name__)
19d2r = num.pi/180.
20r2d = 180./num.pi
21km = 1e3
24class AnalyticalSource(Location):
25 '''
26 Base class for analytical source models.
27 '''
29 name = String.T(
30 optional=True,
31 default='')
33 time = Timestamp.T(
34 default=0.,
35 help='Source origin time',
36 optional=True)
38 vr = Float.T(
39 default=0.,
40 help='Rupture velocity [m/s]',
41 optional=True)
43 @property
44 def northing(self):
45 return self.north_shift
47 @property
48 def easting(self):
49 return self.east_shift
52class AnalyticalRectangularSource(AnalyticalSource):
53 '''
54 Rectangular analytical source model.
56 Coordinates on the source plane are with respect to the origin point given
57 by `(lat, lon, east_shift, north_shift, depth)`.
58 '''
60 strike = Float.T(
61 default=0.0,
62 help='Strike direction in [deg], measured clockwise from north')
64 dip = Float.T(
65 default=90.0,
66 help='Dip angle in [deg], measured downward from horizontal')
68 rake = Float.T(
69 default=0.0,
70 help='Rake angle in [deg], measured counter-clockwise from '
71 'right-horizontal in on-plane view')
73 al1 = Float.T(
74 default=0.,
75 help='Left edge source plane coordinate [m]')
77 al2 = Float.T(
78 default=0.,
79 help='Right edge source plane coordinate [m]')
81 aw1 = Float.T(
82 default=0.,
83 help='Lower edge source plane coordinate [m]')
85 aw2 = Float.T(
86 default=0.,
87 help='Upper edge source plane coordinate [m]')
89 slip = Float.T(
90 default=0.,
91 help='Slip on the rectangular source area [m]',
92 optional=True)
94 @property
95 def length(self):
96 return abs(-self.al1 + self.al2)
98 @property
99 def width(self):
100 return abs(-self.aw1 + self.aw2)
102 @property
103 def area(self):
104 return self.width * self.length
107class OkadaSource(AnalyticalRectangularSource):
108 '''
109 Rectangular Okada source model.
110 '''
112 opening = Float.T(
113 default=0.,
114 help='Opening of the plane in [m]',
115 optional=True)
117 poisson__ = Float.T(
118 default=0.25,
119 help='Poisson\'s ratio :math:`\\nu`',
120 optional=True)
122 lamb__ = Float.T(
123 help='First Lame\' s parameter :math:`\\lambda` [Pa]',
124 optional=True)
126 shearmod__ = Float.T(
127 default=32.0e9,
128 help='Shear modulus along the plane :math:`\\mu` [Pa]',
129 optional=True)
131 @property
132 def poisson(self):
133 '''
134 Calculation of Poisson\' s ratio :math:`\\nu` (if not given).
136 The Poisson\' s ratio :math:`\\nu` can be calculated from the Lame
137 parameters :math:`\\lambda` and :math:`\\mu` using :math:`\\nu =
138 \\frac{\\lambda}{2(\\lambda + \\mu)}` (e.g. Mueller 2007).
139 '''
141 if self.poisson__ is not None:
142 return self.poisson__
144 if self.shearmod__ is None or self.lamb__ is None:
145 raise ValueError('Shearmod and lambda are needed')
147 return (self.lamb__) / (2. * (self.lamb__ + self.shearmod__))
149 @poisson.setter
150 def poisson(self, poisson):
151 self.poisson__ = poisson
153 @property
154 def lamb(self):
155 '''
156 Calculation of first Lame\' s parameter (if not given).
158 Poisson\' s ratio :math:`\\nu` and shear modulus :math:`\\mu` must be
159 available.
160 '''
162 if self.lamb__ is not None:
163 return self.lamb__
165 if self.shearmod__ is None or self.poisson__ is None:
166 raise ValueError('Shearmod and poisson ratio are needed')
168 return (
169 2. * self.poisson__ * self.shearmod__) / (1. - 2. * self.poisson__)
171 @lamb.setter
172 def lamb(self, lamb):
173 self.lamb__ = lamb
175 @property
176 def shearmod(self):
177 '''
178 Calculation of shear modulus :math:`\\mu` (if not given).
180 Poisson ratio\' s :math:`\\nu` must be available.
182 .. important ::
184 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`
186 Through :math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to
187 :math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`
189 '''
191 if self.shearmod__ is not None:
192 return self.shearmod__
194 if self.poisson__ is None:
195 raise ValueError('Poisson ratio is needed')
197 return (8. * (1. + self.poisson__)) / (1. - 2. * self.poisson__)
199 @shearmod.setter
200 def shearmod(self, shearmod):
201 self.shearmod__ = shearmod
203 @property
204 def seismic_moment(self):
205 '''
206 Scalar Seismic moment :math:`M_0`.
208 Code copied from Kite. It disregards the opening (as for now).
209 We assume :math:`M_0 = mu A D`
211 .. important ::
213 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`
215 Through :math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to
216 :math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`
218 :return: Seismic moment release
219 :rtype: float
220 '''
222 mu = self.shearmod
224 disl = 0.
225 if self.slip:
226 disl = self.slip
227 if self.opening:
228 disl = (disl**2 + self.opening**2)**.5
230 return mu * self.area * disl
232 @property
233 def moment_magnitude(self):
234 '''
235 Moment magnitude :math:`M_\\mathrm{w}` from Seismic moment.
237 We assume :math:`M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0) - 10.7`
239 :returns: Moment magnitude
240 :rtype: float
241 '''
242 return mt.moment_to_magnitude(self.seismic_moment)
244 def source_patch(self):
245 '''
246 Build source information array for okada_ext.okada input.
248 :return: array of the source data as input for okada_ext.okada
249 :rtype: :py:class:`numpy.ndarray`, ``(1, 9)``
250 '''
251 return num.array([
252 self.northing,
253 self.easting,
254 self.depth,
255 self.strike,
256 self.dip,
257 self.al1,
258 self.al2,
259 self.aw1,
260 self.aw2])
262 def source_disloc(self):
263 '''
264 Build source dislocation for okada_ext.okada input.
266 :return: array of the source dislocation data as input for
267 okada_ext.okada
268 :rtype: :py:class:`numpy.ndarray`, ``(1, 3)``
269 '''
270 return num.array([
271 num.cos(self.rake * d2r) * self.slip,
272 num.sin(self.rake * d2r) * self.slip,
273 self.opening])
275 def discretize(self, nlength, nwidth, *args, **kwargs):
276 '''
277 Discretize fault into rectilinear grid of fault patches.
279 Fault orientation, slip and elastic parameters are passed to the
280 sub-faults unchanged.
282 :param nlength: Number of patches in strike direction
283 :type nlength: int
284 :param nwidth: Number of patches in down-dip direction
285 :type nwidth: int
286 :return: Discrete fault patches
287 :rtype: list of :py:class:`~pyrocko.modelling.okada.OkadaPatch` objects
288 '''
289 assert nlength > 0
290 assert nwidth > 0
292 il = num.repeat(num.arange(nlength), nwidth)
293 iw = num.tile(num.arange(nwidth), nlength)
295 patch_length = self.length / nlength
296 patch_width = self.width / nwidth
298 al1 = -patch_length / 2.
299 al2 = patch_length / 2.
300 aw1 = -patch_width / 2.
301 aw2 = patch_width / 2.
303 source_points = num.zeros((nlength * nwidth, 3))
304 source_points[:, 0] = il * patch_length + patch_length / 2.
305 source_points[:, 1] = iw * patch_width + patch_width / 2.
307 source_points[:, 0] += self.al1
308 source_points[:, 1] -= self.aw2
310 rotmat = num.asarray(
311 mt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.))
313 source_points_rot = num.dot(rotmat.T, source_points.T).T
314 source_points_rot[:, 0] += self.northing
315 source_points_rot[:, 1] += self.easting
316 source_points_rot[:, 2] += self.depth
318 kwargs = {
319 prop: getattr(self, prop) for prop in self.T.propnames
320 if prop not in [
321 'north_shift', 'east_shift', 'depth',
322 'al1', 'al2', 'aw1', 'aw2']}
324 return (
325 [OkadaPatch(
326 parent=self,
327 ix=src_point[0],
328 iy=src_point[1],
329 north_shift=coord[0],
330 east_shift=coord[1],
331 depth=coord[2],
332 al1=al1, al2=al2, aw1=aw1, aw2=aw2, **kwargs)
333 for src_point, coord in zip(source_points, source_points_rot)],
334 source_points)
337class OkadaPatch(OkadaSource):
339 '''
340 Okada source with additional 2D indexes for bookkeeping.
341 '''
343 ix = Int.T(help='Relative index of the patch in x')
344 iy = Int.T(help='Relative index of the patch in y')
346 def __init__(self, parent=None, *args, **kwargs):
347 OkadaSource.__init__(self, *args, **kwargs)
348 self.parent = parent
351def make_okada_coefficient_matrix(
352 source_patches_list,
353 pure_shear=False,
354 rotate_sdn=True,
355 nthreads=1, variant='normal'):
357 '''
358 Build coefficient matrix for given fault patches.
360 The boundary element method (BEM) for a discretized fault and the
361 determination of the slip distribution from stress drop is based on the
362 relation :math:`stress = coefmat \\cdot displ`. Here the coefficient matrix
363 is built, based on the displacements from Okada's solution and their
364 partial derivatives.
366 :param source_patches_list: Source patches, to be used in BEM.
367 :type source_patches_list: list of
368 :py:class:`~pyrocko.modelling.okada.OkadaSource` objects
369 :param pure_shear: If ``True``, only shear forces are taken into account.
370 :type pure_shear: optional, bool
371 :param rotate_sdn: If ``True``, rotate to strike, dip, normal.
372 :type rotate_sdn: optional, bool
373 :param nthreads: Number of threads.
374 :type nthreads: optional, int
376 :return: Coefficient matrix for all source combinations.
377 :rtype: :py:class:`numpy.ndarray`,
378 ``(len(source_patches_list) * 3, len(source_patches_list) * 3)``
379 '''
381 if variant == 'slow':
382 return _make_okada_coefficient_matrix_slow(
383 source_patches_list, pure_shear, rotate_sdn, nthreads)
385 source_patches = num.array([
386 src.source_patch() for src in source_patches_list])
387 receiver_coords = source_patches[:, :3].copy()
389 npoints = len(source_patches_list)
391 if pure_shear:
392 n_eq = 2
393 else:
394 n_eq = 3
396 coefmat = num.zeros((npoints * 3, npoints * 3))
398 lambda_mean = num.mean([src.lamb for src in source_patches_list])
399 mu_mean = num.mean([src.shearmod for src in source_patches_list])
401 unit_disl = 1.
402 disl_cases = {
403 'strikeslip': {
404 'slip': unit_disl,
405 'opening': 0.,
406 'rake': 0.},
407 'dipslip': {
408 'slip': unit_disl,
409 'opening': 0.,
410 'rake': 90.},
411 'tensileslip': {
412 'slip': 0.,
413 'opening': unit_disl,
414 'rake': 0.}
415 }
417 diag_ind = [0, 4, 8]
418 kron = num.zeros(9)
419 kron[diag_ind] = 1.
421 if variant == 'normal':
422 kron = kron[num.newaxis, num.newaxis, :]
423 else:
424 kron = kron[num.newaxis, :]
426 for idisl, case_type in enumerate([
427 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
428 case = disl_cases[case_type]
429 source_disl = num.array([
430 case['slip'] * num.cos(case['rake'] * d2r),
431 case['slip'] * num.sin(case['rake'] * d2r),
432 case['opening']])
434 if variant == 'normal':
435 results = okada_ext.okada(
436 source_patches,
437 num.tile(source_disl, npoints).reshape(-1, 3),
438 receiver_coords,
439 lambda_mean,
440 mu_mean,
441 nthreads=nthreads,
442 rotate_sdn=int(rotate_sdn),
443 stack_sources=int(variant != 'normal'))
445 eps = 0.5 * (
446 results[:, :, 3:] +
447 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
449 dilatation \
450 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
452 stress_sdn = kron*lambda_mean*dilatation + 2.*mu_mean*eps
453 coefmat[:, idisl::3] = stress_sdn[:, :, (2, 5, 8)]\
454 .reshape(-1, npoints*3).T
455 else:
456 for isrc, source in enumerate(source_patches):
457 results = okada_ext.okada(
458 source[num.newaxis, :],
459 source_disl[num.newaxis, :],
460 receiver_coords,
461 lambda_mean,
462 mu_mean,
463 nthreads=nthreads,
464 rotate_sdn=int(rotate_sdn))
466 eps = 0.5 * (
467 results[:, 3:] +
468 results[:, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
470 dilatation \
471 = num.sum(eps[:, diag_ind], axis=1)[:, num.newaxis]
472 stress_sdn \
473 = kron * lambda_mean * dilatation+2. * mu_mean * eps
475 coefmat[:, isrc*3 + idisl] \
476 = stress_sdn[:, (2, 5, 8)].ravel()
478 if pure_shear:
479 coefmat[2::3, :] = 0.
481 return -coefmat / unit_disl
484def _make_okada_coefficient_matrix_slow(
485 source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1):
487 source_patches = num.array([
488 src.source_patch() for src in source_patches_list])
489 receiver_coords = source_patches[:, :3].copy()
491 npoints = len(source_patches_list)
493 if pure_shear:
494 n_eq = 2
495 else:
496 n_eq = 3
498 coefmat = num.zeros((npoints * 3, npoints * 3))
500 def ned2sdn_rotmat(strike, dip):
501 rotmat = mt.euler_to_matrix(
502 (dip + 180.) * d2r, strike * d2r, 0.).A
503 return rotmat
505 lambda_mean = num.mean([src.lamb for src in source_patches_list])
506 shearmod_mean = num.mean([src.shearmod for src in source_patches_list])
508 unit_disl = 1.
509 disl_cases = {
510 'strikeslip': {
511 'slip': unit_disl,
512 'opening': 0.,
513 'rake': 0.},
514 'dipslip': {
515 'slip': unit_disl,
516 'opening': 0.,
517 'rake': 90.},
518 'tensileslip': {
519 'slip': 0.,
520 'opening': unit_disl,
521 'rake': 0.}
522 }
523 for idisl, case_type in enumerate([
524 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
525 case = disl_cases[case_type]
526 source_disl = num.array([
527 case['slip'] * num.cos(case['rake'] * d2r),
528 case['slip'] * num.sin(case['rake'] * d2r),
529 case['opening']])
531 for isource, source in enumerate(source_patches):
532 results = okada_ext.okada(
533 source[num.newaxis, :].copy(),
534 source_disl[num.newaxis, :].copy(),
535 receiver_coords,
536 lambda_mean,
537 shearmod_mean,
538 nthreads=nthreads,
539 rotate_sdn=int(rotate_sdn))
541 for irec in range(receiver_coords.shape[0]):
542 eps = num.zeros((3, 3))
543 for m in range(3):
544 for n in range(3):
545 eps[m, n] = 0.5 * (
546 results[irec][m * 3 + n + 3] +
547 results[irec][n * 3 + m + 3])
549 stress = num.zeros((3, 3))
550 dilatation = num.sum([eps[i, i] for i in range(3)])
552 for m, n in zip([0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]):
553 if m == n:
554 stress[m, n] = \
555 lambda_mean * \
556 dilatation + \
557 2. * shearmod_mean * \
558 eps[m, n]
560 else:
561 stress[m, n] = \
562 2. * shearmod_mean * \
563 eps[m, n]
564 stress[n, m] = stress[m, n]
566 normal = num.array([0., 0., -1.])
567 for isig in range(3):
568 tension = num.sum(stress[isig, :] * normal)
569 coefmat[irec * n_eq + isig, isource * n_eq + idisl] = \
570 tension / unit_disl
572 return coefmat
575def invert_fault_dislocations_bem(
576 stress_field,
577 coef_mat=None,
578 source_list=None,
579 pure_shear=False,
580 epsilon=None,
581 nthreads=1,
582 **kwargs):
583 '''
584 BEM least squares inversion to get fault dislocations given stress field.
586 Follows least squares inversion approach by Menke (1989) to calculate
587 dislocations on a fault with several segments from a given stress field.
588 The coefficient matrix connecting stresses and displacements of the fault
589 patches can either be specified by the user (``coef_mat``) or it is
590 calculated using the solution of Okada (1992) for a rectangular fault in a
591 homogeneous half space (``source_list``).
593 :param stress_field: Stress change [Pa] for each source patch (as
594 ``stress_field[isource, icomponent]`` where isource indexes the source
595 patch and ``icomponent`` indexes component, ordered (strike, dip,
596 tensile).
597 :type stress_field: :py:class:`numpy.ndarray`, shape ``(nsources, 3)``
599 :param coef_mat: Coefficient matrix connecting source patch
600 dislocations and the stress field.
601 :type coef_mat: optional, :py:class:`numpy.ndarray`, shape
602 ``(len(source_list) * 3, len(source_list) * 3)``
604 :param source_list: list of all source patches to be used for BEM.
605 :type source_list: optional, list of
606 :py:class:`~pyrocko.modelling.okada.OkadaSource` objects
608 :param epsilon: If given, values in ``coef_mat`` smaller than ``epsilon``
609 are set to zero.
610 :type epsilon: optional, float
612 :param nthreads: Number of threads allowed.
613 :type nthreads: int
615 :return: Inverted displacements as ``displacements[isource, icomponent]``
616 where isource indexes the source patch and ``icomponent`` indexes
617 component, ordered (strike, dip, tensile).
618 :rtype: :py:class:`numpy.ndarray`, ``(nsources, 3)``
619 '''
621 if source_list is not None and coef_mat is None:
622 coef_mat = make_okada_coefficient_matrix(
623 source_list, pure_shear=pure_shear, nthreads=nthreads,
624 **kwargs)
626 if epsilon is not None:
627 coef_mat[coef_mat < epsilon] = 0.
629 idx = num.arange(0, coef_mat.shape[0])
630 if pure_shear:
631 idx = idx[idx % 3 != 2]
633 coef_mat_in = coef_mat[idx, :][:, idx]
634 disloc_est = num.zeros(coef_mat.shape[0])
636 if stress_field.ndim == 2:
637 stress_field = stress_field.ravel()
639 threadpool_limits = get_threadpool_limits()
641 with threadpool_limits(limits=nthreads, user_api='blas'):
642 try:
643 disloc_est[idx] = num.linalg.multi_dot([
644 num.linalg.inv(num.dot(coef_mat_in.T, coef_mat_in)),
645 coef_mat_in.T,
646 stress_field[idx]])
647 except num.linalg.LinAlgError as e:
648 logger.warning('Linear inversion failed!')
649 logger.warning(
650 'coef_mat: %s\nstress_field: %s',
651 coef_mat_in, stress_field[idx])
652 raise e
653 return disloc_est.reshape(-1, 3)
656__all__ = [
657 'AnalyticalSource',
658 'AnalyticalRectangularSource',
659 'OkadaSource',
660 'OkadaPatch',
661 'make_okada_coefficient_matrix',
662 'invert_fault_dislocations_bem']