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 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 First Lame\' s parameter :math:`\\lambda` (if not given).
158 Poisson\' s ratio :math:`\\nu` and shear modulus :math:`\\mu` must be
159 available to calculate the first Lame\' s parameter :math:`\\lambda`.
161 .. important ::
163 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`.
165 Through :math:`\\nu = \\frac{\\lambda}{2(\\lambda + \\mu)}` this
166 leads to :math:`\\lambda = \\frac{2 \\mu \\nu}{1-2\\nu}`.
168 '''
170 if self.lamb__ is not None:
171 return self.lamb__
173 if self.shearmod__ is None or self.poisson__ is None:
174 raise ValueError('Shearmod and poisson ratio are needed')
176 return (
177 2. * self.poisson__ * self.shearmod__) / (1. - 2. * self.poisson__)
179 @lamb.setter
180 def lamb(self, lamb):
181 self.lamb__ = lamb
183 @property
184 def shearmod(self):
185 '''
186 Shear modulus :math:`\\mu` (if not given).
188 Poisson ratio\' s :math:`\\nu` must be available.
190 .. important ::
192 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`.
194 Through :math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to
195 :math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`.
197 '''
199 if self.shearmod__ is not None:
200 return self.shearmod__
202 if self.poisson__ is None:
203 raise ValueError('Poisson ratio is needed')
205 return (8. * (1. + self.poisson__)) / (1. - 2. * self.poisson__)
207 @shearmod.setter
208 def shearmod(self, shearmod):
209 self.shearmod__ = shearmod
211 @property
212 def seismic_moment(self):
213 '''
214 Scalar Seismic moment :math:`M_0`.
216 Code copied from Kite. It disregards the opening (as for now).
217 We assume :math:`M_0 = mu A D`.
219 .. important ::
221 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`.
223 Through :math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to
224 :math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`.
226 :return:
227 Seismic moment release.
228 :rtype:
229 float
230 '''
232 mu = self.shearmod
234 disl = 0.
235 if self.slip:
236 disl = self.slip
237 if self.opening:
238 disl = (disl**2 + self.opening**2)**.5
240 return mu * self.area * disl
242 @property
243 def moment_magnitude(self):
244 '''
245 Moment magnitude :math:`M_\\mathrm{w}` from seismic moment.
247 We assume :math:`M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0) - 10.7`.
249 :returns:
250 Moment magnitude.
251 :rtype:
252 float
253 '''
254 return mt.moment_to_magnitude(self.seismic_moment)
256 def source_patch(self):
257 '''
258 Get source location and geometry array for okada_ext.okada input.
260 The values are defined according to Okada (1992).
262 :return:
263 Source data as input for okada_ext.okada. The order is
264 northing [m], easting [m], depth [m], strike [deg], dip [deg],
265 al1 [m], al2 [m], aw1 [m], aw2 [m].
266 :rtype:
267 :py:class:`~numpy.ndarray`: ``(9, )``
268 '''
269 return num.array([
270 self.northing,
271 self.easting,
272 self.depth,
273 self.strike,
274 self.dip,
275 self.al1,
276 self.al2,
277 self.aw1,
278 self.aw2])
280 def source_disloc(self):
281 '''
282 Get source dislocation array for okada_ext.okada input.
284 The given slip is splitted into a strike and an updip part based on the
285 source rake.
287 :return:
288 Source dislocation data as input for okada_ext.okada. The order is
289 dislocation in strike [m], dislocation updip [m], opening [m].
290 :rtype:
291 :py:class:`~numpy.ndarray`: ``(3, )``
292 '''
293 return num.array([
294 num.cos(self.rake * d2r) * self.slip,
295 num.sin(self.rake * d2r) * self.slip,
296 self.opening])
298 def discretize(self, nlength, nwidth, *args, **kwargs):
299 '''
300 Discretize fault into rectilinear grid of fault patches.
302 Fault orientation, slip and elastic parameters are passed to the
303 sub-faults unchanged.
305 :param nlength:
306 Number of patches in strike direction.
307 :type nlength:
308 int
310 :param nwidth:
311 Number of patches in down-dip direction.
312 :type nwidth:
313 int
315 :return:
316 Discrete fault patches.
317 :rtype:
318 list of :py:class:`~pyrocko.modelling.okada.OkadaPatch`
319 '''
320 assert nlength > 0
321 assert nwidth > 0
323 il = num.repeat(num.arange(nlength), nwidth)
324 iw = num.tile(num.arange(nwidth), nlength)
326 patch_length = self.length / nlength
327 patch_width = self.width / nwidth
329 al1 = -patch_length / 2.
330 al2 = patch_length / 2.
331 aw1 = -patch_width / 2.
332 aw2 = patch_width / 2.
334 source_points = num.zeros((nlength * nwidth, 3))
335 source_points[:, 0] = il * patch_length + patch_length / 2.
336 source_points[:, 1] = iw * patch_width + patch_width / 2.
338 source_points[:, 0] += self.al1
339 source_points[:, 1] -= self.aw2
341 rotmat = num.asarray(
342 mt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.))
344 source_points_rot = num.dot(rotmat.T, source_points.T).T
345 source_points_rot[:, 0] += self.northing
346 source_points_rot[:, 1] += self.easting
347 source_points_rot[:, 2] += self.depth
349 kwargs = {
350 prop: getattr(self, prop) for prop in self.T.propnames
351 if prop not in [
352 'north_shift', 'east_shift', 'depth',
353 'al1', 'al2', 'aw1', 'aw2']}
355 return (
356 [OkadaPatch(
357 parent=self,
358 ix=src_point[0],
359 iy=src_point[1],
360 north_shift=coord[0],
361 east_shift=coord[1],
362 depth=coord[2],
363 al1=al1, al2=al2, aw1=aw1, aw2=aw2, **kwargs)
364 for src_point, coord in zip(source_points, source_points_rot)],
365 source_points)
368class OkadaPatch(OkadaSource):
370 '''
371 Okada source with additional 2D indexes for bookkeeping.
372 '''
374 ix = Int.T(help='Relative index of the patch in x')
375 iy = Int.T(help='Relative index of the patch in y')
377 def __init__(self, parent=None, *args, **kwargs):
378 OkadaSource.__init__(self, *args, **kwargs)
379 self.parent = parent
382def make_okada_coefficient_matrix(
383 source_patches_list,
384 pure_shear=False,
385 rotate_sdn=True,
386 nthreads=1, variant='normal'):
388 '''
389 Build coefficient matrix for given fault patches.
391 The boundary element method (BEM) for a discretized fault and the
392 determination of the slip distribution :math:`\\Delta u` from stress drop
393 :math:`\\Delta \\sigma` is based on
394 :math:`\\Delta \\sigma = \\mathbf{C} \\cdot \\Delta u`. Here the
395 coefficient matrix :math:`\\mathbf{C}` is built, based on the displacements
396 from Okada's solution (Okada, 1992) and their partial derivatives.
398 :param source_patches_list:
399 Source patches, to be used in BEM.
400 :type source_patches_list:
401 list of :py:class:`~pyrocko.modelling.okada.OkadaSource`.
403 :param pure_shear:
404 If ``True``, only shear forces are taken into account.
405 :type pure_shear:
406 optional, bool
408 :param rotate_sdn:
409 If ``True``, rotate to strike, dip, normal.
410 :type rotate_sdn:
411 optional, bool
413 :param nthreads:
414 Number of threads.
415 :type nthreads:
416 optional, int
418 :return:
419 Coefficient matrix for all source combinations.
420 :rtype:
421 :py:class:`~numpy.ndarray`:
422 ``(len(source_patches_list) * 3, len(source_patches_list) * 3)``
423 '''
425 if variant == 'slow':
426 return _make_okada_coefficient_matrix_slow(
427 source_patches_list, pure_shear, rotate_sdn, nthreads)
429 source_patches = num.array([
430 src.source_patch() for src in source_patches_list])
431 receiver_coords = source_patches[:, :3].copy()
433 npoints = len(source_patches_list)
435 if pure_shear:
436 n_eq = 2
437 else:
438 n_eq = 3
440 coefmat = num.zeros((npoints * 3, npoints * 3))
442 lambda_mean = num.mean([src.lamb for src in source_patches_list])
443 mu_mean = num.mean([src.shearmod for src in source_patches_list])
445 unit_disl = 1.
446 disl_cases = {
447 'strikeslip': {
448 'slip': unit_disl,
449 'opening': 0.,
450 'rake': 0.},
451 'dipslip': {
452 'slip': unit_disl,
453 'opening': 0.,
454 'rake': 90.},
455 'tensileslip': {
456 'slip': 0.,
457 'opening': unit_disl,
458 'rake': 0.}
459 }
461 diag_ind = [0, 4, 8]
462 kron = num.zeros(9)
463 kron[diag_ind] = 1.
465 if variant == 'normal':
466 kron = kron[num.newaxis, num.newaxis, :]
467 else:
468 kron = kron[num.newaxis, :]
470 for idisl, case_type in enumerate([
471 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
472 case = disl_cases[case_type]
473 source_disl = num.array([
474 case['slip'] * num.cos(case['rake'] * d2r),
475 case['slip'] * num.sin(case['rake'] * d2r),
476 case['opening']])
478 if variant == 'normal':
479 results = okada_ext.okada(
480 source_patches,
481 num.tile(source_disl, npoints).reshape(-1, 3),
482 receiver_coords,
483 lambda_mean,
484 mu_mean,
485 nthreads=nthreads,
486 rotate_sdn=int(rotate_sdn),
487 stack_sources=int(variant != 'normal'))
489 eps = 0.5 * (
490 results[:, :, 3:] +
491 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
493 dilatation \
494 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
496 stress_sdn = kron*lambda_mean*dilatation + 2.*mu_mean*eps
497 coefmat[:, idisl::3] = stress_sdn[:, :, (2, 5, 8)]\
498 .reshape(-1, npoints*3).T
499 else:
500 for isrc, source in enumerate(source_patches):
501 results = okada_ext.okada(
502 source[num.newaxis, :],
503 source_disl[num.newaxis, :],
504 receiver_coords,
505 lambda_mean,
506 mu_mean,
507 nthreads=nthreads,
508 rotate_sdn=int(rotate_sdn))
510 eps = 0.5 * (
511 results[:, 3:] +
512 results[:, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
514 dilatation \
515 = num.sum(eps[:, diag_ind], axis=1)[:, num.newaxis]
516 stress_sdn \
517 = kron * lambda_mean * dilatation+2. * mu_mean * eps
519 coefmat[:, isrc*3 + idisl] \
520 = stress_sdn[:, (2, 5, 8)].ravel()
522 if pure_shear:
523 coefmat[2::3, :] = 0.
525 return -coefmat / unit_disl
528def _make_okada_coefficient_matrix_slow(
529 source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1):
531 source_patches = num.array([
532 src.source_patch() for src in source_patches_list])
533 receiver_coords = source_patches[:, :3].copy()
535 npoints = len(source_patches_list)
537 if pure_shear:
538 n_eq = 2
539 else:
540 n_eq = 3
542 coefmat = num.zeros((npoints * 3, npoints * 3))
544 def ned2sdn_rotmat(strike, dip):
545 rotmat = mt.euler_to_matrix(
546 (dip + 180.) * d2r, strike * d2r, 0.).A
547 return rotmat
549 lambda_mean = num.mean([src.lamb for src in source_patches_list])
550 shearmod_mean = num.mean([src.shearmod for src in source_patches_list])
552 unit_disl = 1.
553 disl_cases = {
554 'strikeslip': {
555 'slip': unit_disl,
556 'opening': 0.,
557 'rake': 0.},
558 'dipslip': {
559 'slip': unit_disl,
560 'opening': 0.,
561 'rake': 90.},
562 'tensileslip': {
563 'slip': 0.,
564 'opening': unit_disl,
565 'rake': 0.}
566 }
567 for idisl, case_type in enumerate([
568 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
569 case = disl_cases[case_type]
570 source_disl = num.array([
571 case['slip'] * num.cos(case['rake'] * d2r),
572 case['slip'] * num.sin(case['rake'] * d2r),
573 case['opening']])
575 for isource, source in enumerate(source_patches):
576 results = okada_ext.okada(
577 source[num.newaxis, :].copy(),
578 source_disl[num.newaxis, :].copy(),
579 receiver_coords,
580 lambda_mean,
581 shearmod_mean,
582 nthreads=nthreads,
583 rotate_sdn=int(rotate_sdn))
585 for irec in range(receiver_coords.shape[0]):
586 eps = num.zeros((3, 3))
587 for m in range(3):
588 for n in range(3):
589 eps[m, n] = 0.5 * (
590 results[irec][m * 3 + n + 3] +
591 results[irec][n * 3 + m + 3])
593 stress = num.zeros((3, 3))
594 dilatation = num.sum([eps[i, i] for i in range(3)])
596 for m, n in zip([0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]):
597 if m == n:
598 stress[m, n] = \
599 lambda_mean * \
600 dilatation + \
601 2. * shearmod_mean * \
602 eps[m, n]
604 else:
605 stress[m, n] = \
606 2. * shearmod_mean * \
607 eps[m, n]
608 stress[n, m] = stress[m, n]
610 normal = num.array([0., 0., -1.])
611 for isig in range(3):
612 tension = num.sum(stress[isig, :] * normal)
613 coefmat[irec * n_eq + isig, isource * n_eq + idisl] = \
614 tension / unit_disl
616 return coefmat
619def invert_fault_dislocations_bem(
620 stress_field,
621 coef_mat=None,
622 source_list=None,
623 pure_shear=False,
624 epsilon=None,
625 nthreads=1,
626 **kwargs):
627 '''
628 BEM least squares inversion to get fault dislocations given stress field.
630 Follows least squares inversion approach by Menke (1989) to calculate
631 dislocations on a fault with several segments from a given stress field.
632 The coefficient matrix connecting stresses and displacements of the fault
633 patches can either be specified by the user (``coef_mat``) or it is
634 calculated using the solution of Okada (1992) for a rectangular fault in a
635 homogeneous half space (``source_list``).
637 :param stress_field:
638 Stress change [Pa] for each source patch (as
639 ``stress_field[isource, icomponent]`` where isource indexes the source
640 patch and ``icomponent`` indexes component, ordered (strike, dip,
641 tensile).
642 :type stress_field:
643 :py:class:`~numpy.ndarray`: ``(nsources, 3)``
645 :param coef_mat:
646 Coefficient matrix connecting source patch dislocations and the stress
647 field.
648 :type coef_mat:
649 optional, :py:class:`~numpy.ndarray`:
650 ``(len(source_list) * 3, len(source_list) * 3)``
652 :param source_list:
653 Source patches to be used for BEM.
654 :type source_list:
655 optional, list of
656 :py:class:`~pyrocko.modelling.okada.OkadaSource`
658 :param epsilon:
659 If given, values in ``coef_mat`` smaller than ``epsilon`` are set to
660 zero.
661 :type epsilon:
662 optional, float
664 :param nthreads:
665 Number of threads allowed.
666 :type nthreads:
667 int
669 :return:
670 Inverted displacements as ``displacements[isource, icomponent]``
671 where isource indexes the source patch and ``icomponent`` indexes
672 component, ordered (strike, dip, tensile).
673 :rtype:
674 :py:class:`~numpy.ndarray`: ``(nsources, 3)``
675 '''
677 if source_list is not None and coef_mat is None:
678 coef_mat = make_okada_coefficient_matrix(
679 source_list, pure_shear=pure_shear, nthreads=nthreads,
680 **kwargs)
682 if epsilon is not None:
683 coef_mat[coef_mat < epsilon] = 0.
685 idx = num.arange(0, coef_mat.shape[0])
686 if pure_shear:
687 idx = idx[idx % 3 != 2]
689 coef_mat_in = coef_mat[idx, :][:, idx]
690 disloc_est = num.zeros(coef_mat.shape[0])
692 if stress_field.ndim == 2:
693 stress_field = stress_field.ravel()
695 threadpool_limits = get_threadpool_limits()
697 with threadpool_limits(limits=nthreads, user_api='blas'):
698 try:
699 disloc_est[idx] = num.linalg.multi_dot([
700 num.linalg.inv(num.dot(coef_mat_in.T, coef_mat_in)),
701 coef_mat_in.T,
702 stress_field[idx]])
703 except num.linalg.LinAlgError as e:
704 logger.warning('Linear inversion failed!')
705 logger.warning(
706 'coef_mat: %s\nstress_field: %s',
707 coef_mat_in, stress_field[idx])
708 raise e
709 return disloc_est.reshape(-1, 3)
712__all__ = [
713 'AnalyticalSource',
714 'AnalyticalRectangularSource',
715 'OkadaSource',
716 'OkadaPatch',
717 'make_okada_coefficient_matrix',
718 'invert_fault_dislocations_bem']