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 ratio :math:`\\nu`. '
120 'The Poisson ratio :math:`\\nu`. If set to ``None``, calculated '
121 'from the Lame\' parameters :math:`\\lambda` and :math:`\\mu` '
122 'using :math:`\\nu = \\frac{\\lambda}{2(\\lambda + \\mu)}` (e.g. '
123 'Mueller 2007).',
124 optional=True)
126 lamb__ = Float.T(
127 help='First Lame parameter :math:`\\lambda` [Pa]. '
128 'If set to ``None``, it is computed from Poisson ratio '
129 ':math:`\\nu` and shear modulus :math:`\\mu`. **Important:** We '
130 'assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`. '
131 'Through :math:`\\nu = \\frac{\\lambda}{2(\\lambda + \\mu)}` '
132 'this leads to :math:`\\lambda = \\frac{2 \\mu \\nu}{1-2\\nu}`.',
133 optional=True)
135 shearmod__ = Float.T(
136 default=32.0e9,
137 help='Shear modulus :math:`\\mu` [Pa]. '
138 'If set to ``None``, it is computed from poisson ratio. '
139 '**Important:** We assume a perfect elastic solid with '
140 ':math:`K=\\frac{5}{3}\\mu`. Through '
141 ':math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to '
142 ':math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`.',
143 optional=True)
145 @property
146 def poisson(self):
147 if self.poisson__ is not None:
148 return self.poisson__
150 if self.shearmod__ is None or self.lamb__ is None:
151 raise ValueError('Shearmod and lambda are needed')
153 return (self.lamb__) / (2. * (self.lamb__ + self.shearmod__))
155 @poisson.setter
156 def poisson(self, poisson):
157 self.poisson__ = poisson
159 @property
160 def lamb(self):
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):
178 if self.shearmod__ is not None:
179 return self.shearmod__
181 if self.poisson__ is None:
182 raise ValueError('Poisson ratio is needed')
184 return (8. * (1. + self.poisson__)) / (1. - 2. * self.poisson__)
186 @shearmod.setter
187 def shearmod(self, shearmod):
188 self.shearmod__ = shearmod
190 @property
191 def seismic_moment(self):
192 '''
193 Scalar Seismic moment :math:`M_0`.
195 Code copied from Kite. It disregards the opening (as for now).
196 We assume :math:`M_0 = mu A D`.
198 .. important ::
200 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`.
202 Through :math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to
203 :math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`.
205 :return:
206 Seismic moment release.
207 :rtype:
208 float
209 '''
211 mu = self.shearmod
213 disl = 0.
214 if self.slip:
215 disl = self.slip
216 if self.opening:
217 disl = (disl**2 + self.opening**2)**.5
219 return mu * self.area * disl
221 @property
222 def moment_magnitude(self):
223 '''
224 Moment magnitude :math:`M_\\mathrm{w}` from seismic moment.
226 We assume :math:`M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0) - 10.7`.
228 :returns:
229 Moment magnitude.
230 :rtype:
231 float
232 '''
233 return mt.moment_to_magnitude(self.seismic_moment)
235 def source_patch(self):
236 '''
237 Get source location and geometry array for okada_ext.okada input.
239 The values are defined according to Okada (1992).
241 :return:
242 Source data as input for okada_ext.okada. The order is
243 northing [m], easting [m], depth [m], strike [deg], dip [deg],
244 al1 [m], al2 [m], aw1 [m], aw2 [m].
245 :rtype:
246 :py:class:`~numpy.ndarray`: ``(9, )``
247 '''
248 return num.array([
249 self.northing,
250 self.easting,
251 self.depth,
252 self.strike,
253 self.dip,
254 self.al1,
255 self.al2,
256 self.aw1,
257 self.aw2])
259 def source_disloc(self):
260 '''
261 Get source dislocation array for okada_ext.okada input.
263 The given slip is splitted into a strike and an updip part based on the
264 source rake.
266 :return:
267 Source dislocation data as input for okada_ext.okada. The order is
268 dislocation in strike [m], dislocation updip [m], opening [m].
269 :rtype:
270 :py:class:`~numpy.ndarray`: ``(3, )``
271 '''
272 return num.array([
273 num.cos(self.rake * d2r) * self.slip,
274 num.sin(self.rake * d2r) * self.slip,
275 self.opening])
277 def discretize(self, nlength, nwidth, *args, **kwargs):
278 '''
279 Discretize fault into rectilinear grid of fault patches.
281 Fault orientation, slip and elastic parameters are passed to the
282 sub-faults unchanged.
284 :param nlength:
285 Number of patches in strike direction.
286 :type nlength:
287 int
289 :param nwidth:
290 Number of patches in down-dip direction.
291 :type nwidth:
292 int
294 :return:
295 Discrete fault patches.
296 :rtype:
297 list of :py:class:`~pyrocko.modelling.okada.OkadaPatch`
298 '''
299 assert nlength > 0
300 assert nwidth > 0
302 il = num.repeat(num.arange(nlength), nwidth)
303 iw = num.tile(num.arange(nwidth), nlength)
305 patch_length = self.length / nlength
306 patch_width = self.width / nwidth
308 al1 = -patch_length / 2.
309 al2 = patch_length / 2.
310 aw1 = -patch_width / 2.
311 aw2 = patch_width / 2.
313 source_points = num.zeros((nlength * nwidth, 3))
314 source_points[:, 0] = il * patch_length + patch_length / 2.
315 source_points[:, 1] = iw * patch_width + patch_width / 2.
317 source_points[:, 0] += self.al1
318 source_points[:, 1] -= self.aw2
320 rotmat = mt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.)
322 source_points_rot = num.dot(rotmat.T, source_points.T).T
323 source_points_rot[:, 0] += self.northing
324 source_points_rot[:, 1] += self.easting
325 source_points_rot[:, 2] += self.depth
327 kwargs = {
328 prop: getattr(self, prop) for prop in self.T.propnames
329 if prop not in [
330 'north_shift', 'east_shift', 'depth',
331 'al1', 'al2', 'aw1', 'aw2']}
333 return (
334 [OkadaPatch(
335 parent=self,
336 ix=src_point[0],
337 iy=src_point[1],
338 north_shift=coord[0],
339 east_shift=coord[1],
340 depth=coord[2],
341 al1=al1, al2=al2, aw1=aw1, aw2=aw2, **kwargs)
342 for src_point, coord in zip(source_points, source_points_rot)],
343 source_points)
346class OkadaPatch(OkadaSource):
348 '''
349 Okada source with additional 2D indexes for bookkeeping.
350 '''
352 ix = Int.T(help='Relative index of the patch in x')
353 iy = Int.T(help='Relative index of the patch in y')
355 def __init__(self, parent=None, *args, **kwargs):
356 OkadaSource.__init__(self, *args, **kwargs)
357 self.parent = parent
360def make_okada_coefficient_matrix(
361 source_patches_list,
362 pure_shear=False,
363 rotate_sdn=True,
364 nthreads=1, variant='normal'):
366 '''
367 Build coefficient matrix for given fault patches.
369 The boundary element method (BEM) for a discretized fault and the
370 determination of the slip distribution :math:`\\Delta u` from stress drop
371 :math:`\\Delta \\sigma` is based on
372 :math:`\\Delta \\sigma = \\mathbf{C} \\cdot \\Delta u`. Here the
373 coefficient matrix :math:`\\mathbf{C}` is built, based on the displacements
374 from Okada's solution (Okada, 1992) and their partial derivatives.
376 :param source_patches_list:
377 Source patches, to be used in BEM.
378 :type source_patches_list:
379 list of :py:class:`~pyrocko.modelling.okada.OkadaSource`.
381 :param pure_shear:
382 If ``True``, only shear forces are taken into account.
383 :type pure_shear:
384 optional, bool
386 :param rotate_sdn:
387 If ``True``, rotate to strike, dip, normal.
388 :type rotate_sdn:
389 optional, bool
391 :param nthreads:
392 Number of threads.
393 :type nthreads:
394 optional, int
396 :return:
397 Coefficient matrix for all source combinations.
398 :rtype:
399 :py:class:`~numpy.ndarray`:
400 ``(len(source_patches_list) * 3, len(source_patches_list) * 3)``
401 '''
403 if variant == 'slow':
404 return _make_okada_coefficient_matrix_slow(
405 source_patches_list, pure_shear, rotate_sdn, nthreads)
407 source_patches = num.array([
408 src.source_patch() for src in source_patches_list])
409 receiver_coords = source_patches[:, :3].copy()
411 npoints = len(source_patches_list)
413 if pure_shear:
414 n_eq = 2
415 else:
416 n_eq = 3
418 coefmat = num.zeros((npoints * 3, npoints * 3))
420 lambda_mean = num.mean([src.lamb for src in source_patches_list])
421 mu_mean = num.mean([src.shearmod for src in source_patches_list])
423 unit_disl = 1.
424 disl_cases = {
425 'strikeslip': {
426 'slip': unit_disl,
427 'opening': 0.,
428 'rake': 0.},
429 'dipslip': {
430 'slip': unit_disl,
431 'opening': 0.,
432 'rake': 90.},
433 'tensileslip': {
434 'slip': 0.,
435 'opening': unit_disl,
436 'rake': 0.}
437 }
439 diag_ind = [0, 4, 8]
440 kron = num.zeros(9)
441 kron[diag_ind] = 1.
443 if variant == 'normal':
444 kron = kron[num.newaxis, num.newaxis, :]
445 else:
446 kron = kron[num.newaxis, :]
448 for idisl, case_type in enumerate([
449 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
450 case = disl_cases[case_type]
451 source_disl = num.array([
452 case['slip'] * num.cos(case['rake'] * d2r),
453 case['slip'] * num.sin(case['rake'] * d2r),
454 case['opening']])
456 if variant == 'normal':
457 results = okada_ext.okada(
458 source_patches,
459 num.tile(source_disl, npoints).reshape(-1, 3),
460 receiver_coords,
461 lambda_mean,
462 mu_mean,
463 nthreads=nthreads,
464 rotate_sdn=int(rotate_sdn),
465 stack_sources=int(variant != 'normal'))
467 eps = 0.5 * (
468 results[:, :, 3:] +
469 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
471 dilatation \
472 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
474 stress_sdn = kron*lambda_mean*dilatation + 2.*mu_mean*eps
475 coefmat[:, idisl::3] = stress_sdn[:, :, (2, 5, 8)]\
476 .reshape(-1, npoints*3).T
477 else:
478 for isrc, source in enumerate(source_patches):
479 results = okada_ext.okada(
480 source.reshape(1, -1),
481 source_disl.reshape(1, -1),
482 receiver_coords,
483 lambda_mean,
484 mu_mean,
485 nthreads=nthreads,
486 rotate_sdn=int(rotate_sdn))
488 eps = 0.5 * (
489 results[:, 3:] +
490 results[:, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
492 dilatation \
493 = num.sum(eps[:, diag_ind], axis=1)[:, num.newaxis]
494 stress_sdn \
495 = kron * lambda_mean * dilatation+2. * mu_mean * eps
497 coefmat[:, isrc*3 + idisl] \
498 = stress_sdn[:, (2, 5, 8)].ravel()
500 if pure_shear:
501 coefmat[2::3, :] = 0.
503 return -coefmat / unit_disl
506def _make_okada_coefficient_matrix_slow(
507 source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1):
509 source_patches = num.array([
510 src.source_patch() for src in source_patches_list])
511 receiver_coords = source_patches[:, :3].copy()
513 npoints = len(source_patches_list)
515 if pure_shear:
516 n_eq = 2
517 else:
518 n_eq = 3
520 coefmat = num.zeros((npoints * 3, npoints * 3))
522 def ned2sdn_rotmat(strike, dip):
523 rotmat = mt.euler_to_matrix(
524 (dip + 180.) * d2r, strike * d2r, 0.)
525 return rotmat
527 lambda_mean = num.mean([src.lamb for src in source_patches_list])
528 shearmod_mean = num.mean([src.shearmod for src in source_patches_list])
530 unit_disl = 1.
531 disl_cases = {
532 'strikeslip': {
533 'slip': unit_disl,
534 'opening': 0.,
535 'rake': 0.},
536 'dipslip': {
537 'slip': unit_disl,
538 'opening': 0.,
539 'rake': 90.},
540 'tensileslip': {
541 'slip': 0.,
542 'opening': unit_disl,
543 'rake': 0.}
544 }
545 for idisl, case_type in enumerate([
546 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
547 case = disl_cases[case_type]
548 source_disl = num.array([
549 case['slip'] * num.cos(case['rake'] * d2r),
550 case['slip'] * num.sin(case['rake'] * d2r),
551 case['opening']])
553 for isource, source in enumerate(source_patches):
554 results = okada_ext.okada(
555 source[num.newaxis, :].copy(),
556 source_disl[num.newaxis, :].copy(),
557 receiver_coords,
558 lambda_mean,
559 shearmod_mean,
560 nthreads=nthreads,
561 rotate_sdn=int(rotate_sdn))
563 for irec in range(receiver_coords.shape[0]):
564 eps = num.zeros((3, 3))
565 for m in range(3):
566 for n in range(3):
567 eps[m, n] = 0.5 * (
568 results[irec][m * 3 + n + 3] +
569 results[irec][n * 3 + m + 3])
571 stress = num.zeros((3, 3))
572 dilatation = num.sum([eps[i, i] for i in range(3)])
574 for m, n in zip([0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]):
575 if m == n:
576 stress[m, n] = \
577 lambda_mean * \
578 dilatation + \
579 2. * shearmod_mean * \
580 eps[m, n]
582 else:
583 stress[m, n] = \
584 2. * shearmod_mean * \
585 eps[m, n]
586 stress[n, m] = stress[m, n]
588 normal = num.array([0., 0., -1.])
589 for isig in range(3):
590 tension = num.sum(stress[isig, :] * normal)
591 coefmat[irec * n_eq + isig, isource * n_eq + idisl] = \
592 tension / unit_disl
594 return coefmat
597def invert_fault_dislocations_bem(
598 stress_field,
599 coef_mat=None,
600 source_list=None,
601 pure_shear=False,
602 epsilon=None,
603 nthreads=1,
604 **kwargs):
605 '''
606 BEM least squares inversion to get fault dislocations given stress field.
608 Follows least squares inversion approach by Menke (1989) to calculate
609 dislocations on a fault with several segments from a given stress field.
610 The coefficient matrix connecting stresses and displacements of the fault
611 patches can either be specified by the user (``coef_mat``) or it is
612 calculated using the solution of Okada (1992) for a rectangular fault in a
613 homogeneous half space (``source_list``).
615 :param stress_field:
616 Stress change [Pa] for each source patch (as
617 ``stress_field[isource, icomponent]`` where isource indexes the source
618 patch and ``icomponent`` indexes component, ordered (strike, dip,
619 tensile).
620 :type stress_field:
621 :py:class:`~numpy.ndarray`: ``(nsources, 3)``
623 :param coef_mat:
624 Coefficient matrix connecting source patch dislocations and the stress
625 field.
626 :type coef_mat:
627 optional, :py:class:`~numpy.ndarray`:
628 ``(len(source_list) * 3, len(source_list) * 3)``
630 :param source_list:
631 Source patches to be used for BEM.
632 :type source_list:
633 optional, list of
634 :py:class:`~pyrocko.modelling.okada.OkadaSource`
636 :param epsilon:
637 If given, values in ``coef_mat`` smaller than ``epsilon`` are set to
638 zero.
639 :type epsilon:
640 optional, float
642 :param nthreads:
643 Number of threads allowed.
644 :type nthreads:
645 int
647 :return:
648 Inverted displacements as ``displacements[isource, icomponent]``
649 where isource indexes the source patch and ``icomponent`` indexes
650 component, ordered (strike, dip, tensile).
651 :rtype:
652 :py:class:`~numpy.ndarray`: ``(nsources, 3)``
653 '''
655 if source_list is not None and coef_mat is None:
656 coef_mat = make_okada_coefficient_matrix(
657 source_list, pure_shear=pure_shear, nthreads=nthreads,
658 **kwargs)
660 if epsilon is not None:
661 coef_mat[coef_mat < epsilon] = 0.
663 idx = num.arange(0, coef_mat.shape[0])
664 if pure_shear:
665 idx = idx[idx % 3 != 2]
667 coef_mat_in = coef_mat[idx, :][:, idx]
668 disloc_est = num.zeros(coef_mat.shape[0])
670 if stress_field.ndim == 2:
671 stress_field = stress_field.ravel()
673 threadpool_limits = get_threadpool_limits()
675 with threadpool_limits(limits=nthreads, user_api='blas'):
676 try:
677 disloc_est[idx] = num.linalg.multi_dot([
678 num.linalg.inv(num.dot(coef_mat_in.T, coef_mat_in)),
679 coef_mat_in.T,
680 stress_field[idx]])
681 except num.linalg.LinAlgError as e:
682 logger.warning('Linear inversion failed!')
683 logger.warning(
684 'coef_mat: %s\nstress_field: %s',
685 coef_mat_in, stress_field[idx])
686 raise e
687 return disloc_est.reshape(-1, 3)
690__all__ = [
691 'AnalyticalSource',
692 'AnalyticalRectangularSource',
693 'OkadaSource',
694 'OkadaPatch',
695 'make_okada_coefficient_matrix',
696 'invert_fault_dislocations_bem']