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:
219 Seismic moment release.
220 :rtype:
221 float
222 '''
224 mu = self.shearmod
226 disl = 0.
227 if self.slip:
228 disl = self.slip
229 if self.opening:
230 disl = (disl**2 + self.opening**2)**.5
232 return mu * self.area * disl
234 @property
235 def moment_magnitude(self):
236 '''
237 Moment magnitude :math:`M_\\mathrm{w}` from seismic moment.
239 We assume :math:`M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0) - 10.7`.
241 :returns:
242 Moment magnitude.
243 :rtype:
244 float
245 '''
246 return mt.moment_to_magnitude(self.seismic_moment)
248 def source_patch(self):
249 '''
250 Build source information array for okada_ext.okada input.
252 :return:
253 Source data as input for okada_ext.okada.
254 :rtype:
255 :py:class:`~numpy.ndarray`: ``(9, )``
256 '''
257 return num.array([
258 self.northing,
259 self.easting,
260 self.depth,
261 self.strike,
262 self.dip,
263 self.al1,
264 self.al2,
265 self.aw1,
266 self.aw2])
268 def source_disloc(self):
269 '''
270 Build source dislocation for okada_ext.okada input.
272 :return:
273 Source dislocation data as input for okada_ext.okada
274 :rtype:
275 :py:class:`~numpy.ndarray`: ``(3, )``
276 '''
277 return num.array([
278 num.cos(self.rake * d2r) * self.slip,
279 num.sin(self.rake * d2r) * self.slip,
280 self.opening])
282 def discretize(self, nlength, nwidth, *args, **kwargs):
283 '''
284 Discretize fault into rectilinear grid of fault patches.
286 Fault orientation, slip and elastic parameters are passed to the
287 sub-faults unchanged.
289 :param nlength:
290 Number of patches in strike direction.
291 :type nlength:
292 int
294 :param nwidth:
295 Number of patches in down-dip direction.
296 :type nwidth:
297 int
299 :return:
300 Discrete fault patches.
301 :rtype:
302 list of :py:class:`~pyrocko.modelling.okada.OkadaPatch`
303 '''
304 assert nlength > 0
305 assert nwidth > 0
307 il = num.repeat(num.arange(nlength), nwidth)
308 iw = num.tile(num.arange(nwidth), nlength)
310 patch_length = self.length / nlength
311 patch_width = self.width / nwidth
313 al1 = -patch_length / 2.
314 al2 = patch_length / 2.
315 aw1 = -patch_width / 2.
316 aw2 = patch_width / 2.
318 source_points = num.zeros((nlength * nwidth, 3))
319 source_points[:, 0] = il * patch_length + patch_length / 2.
320 source_points[:, 1] = iw * patch_width + patch_width / 2.
322 source_points[:, 0] += self.al1
323 source_points[:, 1] -= self.aw2
325 rotmat = num.asarray(
326 mt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.))
328 source_points_rot = num.dot(rotmat.T, source_points.T).T
329 source_points_rot[:, 0] += self.northing
330 source_points_rot[:, 1] += self.easting
331 source_points_rot[:, 2] += self.depth
333 kwargs = {
334 prop: getattr(self, prop) for prop in self.T.propnames
335 if prop not in [
336 'north_shift', 'east_shift', 'depth',
337 'al1', 'al2', 'aw1', 'aw2']}
339 return (
340 [OkadaPatch(
341 parent=self,
342 ix=src_point[0],
343 iy=src_point[1],
344 north_shift=coord[0],
345 east_shift=coord[1],
346 depth=coord[2],
347 al1=al1, al2=al2, aw1=aw1, aw2=aw2, **kwargs)
348 for src_point, coord in zip(source_points, source_points_rot)],
349 source_points)
352class OkadaPatch(OkadaSource):
354 '''
355 Okada source with additional 2D indexes for bookkeeping.
356 '''
358 ix = Int.T(help='Relative index of the patch in x')
359 iy = Int.T(help='Relative index of the patch in y')
361 def __init__(self, parent=None, *args, **kwargs):
362 OkadaSource.__init__(self, *args, **kwargs)
363 self.parent = parent
366def make_okada_coefficient_matrix(
367 source_patches_list,
368 pure_shear=False,
369 rotate_sdn=True,
370 nthreads=1, variant='normal'):
372 '''
373 Build coefficient matrix for given fault patches.
375 The boundary element method (BEM) for a discretized fault and the
376 determination of the slip distribution from stress drop is based on the
377 relation :math:`stress = coefmat \\cdot displ`. Here the coefficient matrix
378 is built, based on the displacements from Okada's solution and their
379 partial derivatives.
381 :param source_patches_list:
382 Source patches, to be used in BEM.
383 :type source_patches_list:
384 list of :py:class:`~pyrocko.modelling.okada.OkadaSource`.
385 :param pure_shear:
386 If ``True``, only shear forces are taken into account.
387 :type pure_shear:
388 optional, bool
389 :param rotate_sdn:
390 If ``True``, rotate to strike, dip, normal.
391 :type rotate_sdn:
392 optional, bool
393 :param nthreads:
394 Number of threads.
395 :type nthreads:
396 optional, int
398 :return:
399 Coefficient matrix for all source combinations.
400 :rtype:
401 :py:class:`~numpy.ndarray`:
402 ``(len(source_patches_list) * 3, len(source_patches_list) * 3)``
403 '''
405 if variant == 'slow':
406 return _make_okada_coefficient_matrix_slow(
407 source_patches_list, pure_shear, rotate_sdn, nthreads)
409 source_patches = num.array([
410 src.source_patch() for src in source_patches_list])
411 receiver_coords = source_patches[:, :3].copy()
413 npoints = len(source_patches_list)
415 if pure_shear:
416 n_eq = 2
417 else:
418 n_eq = 3
420 coefmat = num.zeros((npoints * 3, npoints * 3))
422 lambda_mean = num.mean([src.lamb for src in source_patches_list])
423 mu_mean = num.mean([src.shearmod for src in source_patches_list])
425 unit_disl = 1.
426 disl_cases = {
427 'strikeslip': {
428 'slip': unit_disl,
429 'opening': 0.,
430 'rake': 0.},
431 'dipslip': {
432 'slip': unit_disl,
433 'opening': 0.,
434 'rake': 90.},
435 'tensileslip': {
436 'slip': 0.,
437 'opening': unit_disl,
438 'rake': 0.}
439 }
441 diag_ind = [0, 4, 8]
442 kron = num.zeros(9)
443 kron[diag_ind] = 1.
445 if variant == 'normal':
446 kron = kron[num.newaxis, num.newaxis, :]
447 else:
448 kron = kron[num.newaxis, :]
450 for idisl, case_type in enumerate([
451 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
452 case = disl_cases[case_type]
453 source_disl = num.array([
454 case['slip'] * num.cos(case['rake'] * d2r),
455 case['slip'] * num.sin(case['rake'] * d2r),
456 case['opening']])
458 if variant == 'normal':
459 results = okada_ext.okada(
460 source_patches,
461 num.tile(source_disl, npoints).reshape(-1, 3),
462 receiver_coords,
463 lambda_mean,
464 mu_mean,
465 nthreads=nthreads,
466 rotate_sdn=int(rotate_sdn),
467 stack_sources=int(variant != 'normal'))
469 eps = 0.5 * (
470 results[:, :, 3:] +
471 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
473 dilatation \
474 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
476 stress_sdn = kron*lambda_mean*dilatation + 2.*mu_mean*eps
477 coefmat[:, idisl::3] = stress_sdn[:, :, (2, 5, 8)]\
478 .reshape(-1, npoints*3).T
479 else:
480 for isrc, source in enumerate(source_patches):
481 results = okada_ext.okada(
482 source[num.newaxis, :],
483 source_disl[num.newaxis, :],
484 receiver_coords,
485 lambda_mean,
486 mu_mean,
487 nthreads=nthreads,
488 rotate_sdn=int(rotate_sdn))
490 eps = 0.5 * (
491 results[:, 3:] +
492 results[:, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
494 dilatation \
495 = num.sum(eps[:, diag_ind], axis=1)[:, num.newaxis]
496 stress_sdn \
497 = kron * lambda_mean * dilatation+2. * mu_mean * eps
499 coefmat[:, isrc*3 + idisl] \
500 = stress_sdn[:, (2, 5, 8)].ravel()
502 if pure_shear:
503 coefmat[2::3, :] = 0.
505 return -coefmat / unit_disl
508def _make_okada_coefficient_matrix_slow(
509 source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1):
511 source_patches = num.array([
512 src.source_patch() for src in source_patches_list])
513 receiver_coords = source_patches[:, :3].copy()
515 npoints = len(source_patches_list)
517 if pure_shear:
518 n_eq = 2
519 else:
520 n_eq = 3
522 coefmat = num.zeros((npoints * 3, npoints * 3))
524 def ned2sdn_rotmat(strike, dip):
525 rotmat = mt.euler_to_matrix(
526 (dip + 180.) * d2r, strike * d2r, 0.).A
527 return rotmat
529 lambda_mean = num.mean([src.lamb for src in source_patches_list])
530 shearmod_mean = num.mean([src.shearmod for src in source_patches_list])
532 unit_disl = 1.
533 disl_cases = {
534 'strikeslip': {
535 'slip': unit_disl,
536 'opening': 0.,
537 'rake': 0.},
538 'dipslip': {
539 'slip': unit_disl,
540 'opening': 0.,
541 'rake': 90.},
542 'tensileslip': {
543 'slip': 0.,
544 'opening': unit_disl,
545 'rake': 0.}
546 }
547 for idisl, case_type in enumerate([
548 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
549 case = disl_cases[case_type]
550 source_disl = num.array([
551 case['slip'] * num.cos(case['rake'] * d2r),
552 case['slip'] * num.sin(case['rake'] * d2r),
553 case['opening']])
555 for isource, source in enumerate(source_patches):
556 results = okada_ext.okada(
557 source[num.newaxis, :].copy(),
558 source_disl[num.newaxis, :].copy(),
559 receiver_coords,
560 lambda_mean,
561 shearmod_mean,
562 nthreads=nthreads,
563 rotate_sdn=int(rotate_sdn))
565 for irec in range(receiver_coords.shape[0]):
566 eps = num.zeros((3, 3))
567 for m in range(3):
568 for n in range(3):
569 eps[m, n] = 0.5 * (
570 results[irec][m * 3 + n + 3] +
571 results[irec][n * 3 + m + 3])
573 stress = num.zeros((3, 3))
574 dilatation = num.sum([eps[i, i] for i in range(3)])
576 for m, n in zip([0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]):
577 if m == n:
578 stress[m, n] = \
579 lambda_mean * \
580 dilatation + \
581 2. * shearmod_mean * \
582 eps[m, n]
584 else:
585 stress[m, n] = \
586 2. * shearmod_mean * \
587 eps[m, n]
588 stress[n, m] = stress[m, n]
590 normal = num.array([0., 0., -1.])
591 for isig in range(3):
592 tension = num.sum(stress[isig, :] * normal)
593 coefmat[irec * n_eq + isig, isource * n_eq + idisl] = \
594 tension / unit_disl
596 return coefmat
599def invert_fault_dislocations_bem(
600 stress_field,
601 coef_mat=None,
602 source_list=None,
603 pure_shear=False,
604 epsilon=None,
605 nthreads=1,
606 **kwargs):
607 '''
608 BEM least squares inversion to get fault dislocations given stress field.
610 Follows least squares inversion approach by Menke (1989) to calculate
611 dislocations on a fault with several segments from a given stress field.
612 The coefficient matrix connecting stresses and displacements of the fault
613 patches can either be specified by the user (``coef_mat``) or it is
614 calculated using the solution of Okada (1992) for a rectangular fault in a
615 homogeneous half space (``source_list``).
617 :param stress_field:
618 Stress change [Pa] for each source patch (as
619 ``stress_field[isource, icomponent]`` where isource indexes the source
620 patch and ``icomponent`` indexes component, ordered (strike, dip,
621 tensile).
622 :type stress_field:
623 :py:class:`~numpy.ndarray`: ``(nsources, 3)``
625 :param coef_mat:
626 Coefficient matrix connecting source patch dislocations and the stress
627 field.
628 :type coef_mat:
629 optional, :py:class:`~numpy.ndarray`:
630 ``(len(source_list) * 3, len(source_list) * 3)``
632 :param source_list:
633 Source patches to be used for BEM.
634 :type source_list:
635 optional, list of
636 :py:class:`~pyrocko.modelling.okada.OkadaSource`
638 :param epsilon:
639 If given, values in ``coef_mat`` smaller than ``epsilon`` are set to
640 zero.
641 :type epsilon:
642 optional, float
644 :param nthreads:
645 Number of threads allowed.
646 :type nthreads:
647 int
649 :return:
650 Inverted displacements as ``displacements[isource, icomponent]``
651 where isource indexes the source patch and ``icomponent`` indexes
652 component, ordered (strike, dip, tensile).
653 :rtype:
654 :py:class:`~numpy.ndarray`: ``(nsources, 3)``
655 '''
657 if source_list is not None and coef_mat is None:
658 coef_mat = make_okada_coefficient_matrix(
659 source_list, pure_shear=pure_shear, nthreads=nthreads,
660 **kwargs)
662 if epsilon is not None:
663 coef_mat[coef_mat < epsilon] = 0.
665 idx = num.arange(0, coef_mat.shape[0])
666 if pure_shear:
667 idx = idx[idx % 3 != 2]
669 coef_mat_in = coef_mat[idx, :][:, idx]
670 disloc_est = num.zeros(coef_mat.shape[0])
672 if stress_field.ndim == 2:
673 stress_field = stress_field.ravel()
675 threadpool_limits = get_threadpool_limits()
677 with threadpool_limits(limits=nthreads, user_api='blas'):
678 try:
679 disloc_est[idx] = num.linalg.multi_dot([
680 num.linalg.inv(num.dot(coef_mat_in.T, coef_mat_in)),
681 coef_mat_in.T,
682 stress_field[idx]])
683 except num.linalg.LinAlgError as e:
684 logger.warning('Linear inversion failed!')
685 logger.warning(
686 'coef_mat: %s\nstress_field: %s',
687 coef_mat_in, stress_field[idx])
688 raise e
689 return disloc_est.reshape(-1, 3)
692__all__ = [
693 'AnalyticalSource',
694 'AnalyticalRectangularSource',
695 'OkadaSource',
696 'OkadaPatch',
697 'make_okada_coefficient_matrix',
698 'invert_fault_dislocations_bem']