Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/modelling/okada.py: 80%
214 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-04 09:52 +0000
1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Elastostatic solutions and boundary element modelling for rectangular
8dislocation sources.
9'''
11import numpy as num
12import logging
14from pyrocko import moment_tensor as mt
15from pyrocko.guts import Float, String, Timestamp, Int
16from pyrocko.model import Location
17from pyrocko.modelling import okada_ext
18from pyrocko.util import get_threadpool_limits
20guts_prefix = 'modelling'
22logger = logging.getLogger(__name__)
24d2r = num.pi/180.
25r2d = 180./num.pi
26km = 1e3
29class AnalyticalSource(Location):
30 '''
31 Base class for analytical source models.
32 '''
34 name = String.T(
35 optional=True,
36 default='')
38 time = Timestamp.T(
39 default=0.,
40 help='Source origin time',
41 optional=True)
43 vr = Float.T(
44 default=0.,
45 help='Rupture velocity [m/s]',
46 optional=True)
48 @property
49 def northing(self):
50 return self.north_shift
52 @property
53 def easting(self):
54 return self.east_shift
57class AnalyticalRectangularSource(AnalyticalSource):
58 '''
59 Rectangular analytical source model.
61 Coordinates on the source plane are with respect to the origin point given
62 by `(lat, lon, east_shift, north_shift, depth)`.
63 '''
65 strike = Float.T(
66 default=0.0,
67 help='Strike direction in [deg], measured clockwise from north.')
69 dip = Float.T(
70 default=90.0,
71 help='Dip angle in [deg], measured downward from horizontal.')
73 rake = Float.T(
74 default=0.0,
75 help='Rake angle in [deg], measured counter-clockwise from '
76 'right-horizontal in on-plane view.')
78 al1 = Float.T(
79 default=0.,
80 help='Left edge source plane coordinate [m].')
82 al2 = Float.T(
83 default=0.,
84 help='Right edge source plane coordinate [m].')
86 aw1 = Float.T(
87 default=0.,
88 help='Lower edge source plane coordinate [m].')
90 aw2 = Float.T(
91 default=0.,
92 help='Upper edge source plane coordinate [m].')
94 slip = Float.T(
95 default=0.,
96 help='Slip on the rectangular source area [m].',
97 optional=True)
99 @property
100 def length(self):
101 return abs(-self.al1 + self.al2)
103 @property
104 def width(self):
105 return abs(-self.aw1 + self.aw2)
107 @property
108 def area(self):
109 return self.width * self.length
112class OkadaSource(AnalyticalRectangularSource):
113 '''
114 Rectangular Okada source model.
115 '''
117 opening = Float.T(
118 default=0.,
119 help='Opening of the plane in [m].',
120 optional=True)
122 poisson__ = Float.T(
123 default=0.25,
124 help='Poisson ratio :math:`\\nu`. '
125 'The Poisson ratio :math:`\\nu`. If set to ``None``, calculated '
126 "from the Lame' parameters :math:`\\lambda` and :math:`\\mu` "
127 'using :math:`\\nu = \\frac{\\lambda}{2(\\lambda + \\mu)}` (e.g. '
128 'Mueller 2007).',
129 optional=True)
131 lamb__ = Float.T(
132 help='First Lame parameter :math:`\\lambda` [Pa]. '
133 'If set to ``None``, it is computed from Poisson ratio '
134 ':math:`\\nu` and shear modulus :math:`\\mu`. **Important:** We '
135 'assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`. '
136 'Through :math:`\\nu = \\frac{\\lambda}{2(\\lambda + \\mu)}` '
137 'this leads to :math:`\\lambda = \\frac{2 \\mu \\nu}{1-2\\nu}`.',
138 optional=True)
140 shearmod__ = Float.T(
141 default=32.0e9,
142 help='Shear modulus :math:`\\mu` [Pa]. '
143 'If set to ``None``, it is computed from poisson ratio. '
144 '**Important:** We assume a perfect elastic solid with '
145 ':math:`K=\\frac{5}{3}\\mu`. Through '
146 ':math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to '
147 ':math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`.',
148 optional=True)
150 @property
151 def poisson(self):
152 if self.poisson__ is not None:
153 return self.poisson__
155 if self.shearmod__ is None or self.lamb__ is None:
156 raise ValueError('Shearmod and lambda are needed')
158 return (self.lamb__) / (2. * (self.lamb__ + self.shearmod__))
160 @poisson.setter
161 def poisson(self, poisson):
162 self.poisson__ = poisson
164 @property
165 def lamb(self):
167 if self.lamb__ is not None:
168 return self.lamb__
170 if self.shearmod__ is None or self.poisson__ is None:
171 raise ValueError('Shearmod and poisson ratio are needed')
173 return (
174 2. * self.poisson__ * self.shearmod__) / (1. - 2. * self.poisson__)
176 @lamb.setter
177 def lamb(self, lamb):
178 self.lamb__ = lamb
180 @property
181 def shearmod(self):
183 if self.shearmod__ is not None:
184 return self.shearmod__
186 if self.poisson__ is None:
187 raise ValueError('Poisson ratio is needed')
189 return (8. * (1. + self.poisson__)) / (1. - 2. * self.poisson__)
191 @shearmod.setter
192 def shearmod(self, shearmod):
193 self.shearmod__ = shearmod
195 @property
196 def seismic_moment(self):
197 '''
198 Scalar Seismic moment :math:`M_0`.
200 Code copied from Kite. It disregards the opening (as for now).
201 We assume :math:`M_0 = mu A D`.
203 .. important ::
205 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`.
207 Through :math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to
208 :math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`.
210 :return:
211 Seismic moment release.
212 :rtype:
213 float
214 '''
216 mu = self.shearmod
218 disl = 0.
219 if self.slip:
220 disl = self.slip
221 if self.opening:
222 disl = (disl**2 + self.opening**2)**.5
224 return mu * self.area * disl
226 @property
227 def moment_magnitude(self):
228 '''
229 Moment magnitude :math:`M_\\mathrm{w}` from seismic moment.
231 We assume :math:`M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0) - 10.7`.
233 :returns:
234 Moment magnitude.
235 :rtype:
236 float
237 '''
238 return mt.moment_to_magnitude(self.seismic_moment)
240 def source_patch(self):
241 '''
242 Get source location and geometry array for okada_ext.okada input.
244 The values are defined according to Okada (1992).
246 :return:
247 Source data as input for okada_ext.okada. The order is
248 northing [m], easting [m], depth [m], strike [deg], dip [deg],
249 al1 [m], al2 [m], aw1 [m], aw2 [m].
250 :rtype:
251 :py:class:`~numpy.ndarray`: ``(9, )``
252 '''
253 return num.array([
254 self.northing,
255 self.easting,
256 self.depth,
257 self.strike,
258 self.dip,
259 self.al1,
260 self.al2,
261 self.aw1,
262 self.aw2])
264 def source_disloc(self):
265 '''
266 Get source dislocation array for okada_ext.okada input.
268 The given slip is splitted into a strike and an updip part based on the
269 source rake.
271 :return:
272 Source dislocation data as input for okada_ext.okada. The order is
273 dislocation in strike [m], dislocation updip [m], opening [m].
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 = mt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.)
327 source_points_rot = num.dot(rotmat.T, source_points.T).T
328 source_points_rot[:, 0] += self.northing
329 source_points_rot[:, 1] += self.easting
330 source_points_rot[:, 2] += self.depth
332 kwargs = {
333 prop: getattr(self, prop) for prop in self.T.propnames
334 if prop not in [
335 'north_shift', 'east_shift', 'depth',
336 'al1', 'al2', 'aw1', 'aw2']}
338 return (
339 [OkadaPatch(
340 parent=self,
341 ix=src_point[0],
342 iy=src_point[1],
343 north_shift=coord[0],
344 east_shift=coord[1],
345 depth=coord[2],
346 al1=al1, al2=al2, aw1=aw1, aw2=aw2, **kwargs)
347 for src_point, coord in zip(source_points, source_points_rot)],
348 source_points)
351class OkadaPatch(OkadaSource):
353 '''
354 Okada source with additional 2D indexes for bookkeeping.
355 '''
357 ix = Int.T(help='Relative index of the patch in x')
358 iy = Int.T(help='Relative index of the patch in y')
360 def __init__(self, parent=None, *args, **kwargs):
361 OkadaSource.__init__(self, *args, **kwargs)
362 self.parent = parent
365def make_okada_coefficient_matrix(
366 source_patches_list,
367 pure_shear=False,
368 rotate_sdn=True,
369 nthreads=1, variant='normal'):
371 '''
372 Build coefficient matrix for given fault patches.
374 The boundary element method (BEM) for a discretized fault and the
375 determination of the slip distribution :math:`\\Delta u` from stress drop
376 :math:`\\Delta \\sigma` is based on
377 :math:`\\Delta \\sigma = \\mathbf{C} \\cdot \\Delta u`. Here the
378 coefficient matrix :math:`\\mathbf{C}` is built, based on the displacements
379 from Okada's solution (Okada, 1992) and their 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`.
386 :param pure_shear:
387 If ``True``, only shear forces are taken into account.
388 :type pure_shear:
389 bool
391 :param rotate_sdn:
392 If ``True``, rotate to strike, dip, normal.
393 :type rotate_sdn:
394 bool
396 :param nthreads:
397 Number of threads.
398 :type nthreads:
399 int
401 :return:
402 Coefficient matrix for all source combinations.
403 :rtype:
404 :py:class:`~numpy.ndarray`:
405 ``(len(source_patches_list) * 3, len(source_patches_list) * 3)``
406 '''
408 if variant == 'slow':
409 return _make_okada_coefficient_matrix_slow(
410 source_patches_list, pure_shear, rotate_sdn, nthreads)
412 source_patches = num.array([
413 src.source_patch() for src in source_patches_list])
414 receiver_coords = source_patches[:, :3].copy()
416 npoints = len(source_patches_list)
418 if pure_shear:
419 n_eq = 2
420 else:
421 n_eq = 3
423 coefmat = num.zeros((npoints * 3, npoints * 3))
425 lambda_mean = num.mean([src.lamb for src in source_patches_list])
426 mu_mean = num.mean([src.shearmod for src in source_patches_list])
428 unit_disl = 1.
429 disl_cases = {
430 'strikeslip': {
431 'slip': unit_disl,
432 'opening': 0.,
433 'rake': 0.},
434 'dipslip': {
435 'slip': unit_disl,
436 'opening': 0.,
437 'rake': 90.},
438 'tensileslip': {
439 'slip': 0.,
440 'opening': unit_disl,
441 'rake': 0.}
442 }
444 diag_ind = [0, 4, 8]
445 kron = num.zeros(9)
446 kron[diag_ind] = 1.
448 if variant == 'normal':
449 kron = kron[num.newaxis, num.newaxis, :]
450 else:
451 kron = kron[num.newaxis, :]
453 for idisl, case_type in enumerate([
454 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
455 case = disl_cases[case_type]
456 source_disl = num.array([
457 case['slip'] * num.cos(case['rake'] * d2r),
458 case['slip'] * num.sin(case['rake'] * d2r),
459 case['opening']])
461 if variant == 'normal':
462 results = okada_ext.okada(
463 source_patches,
464 num.tile(source_disl, npoints).reshape(-1, 3),
465 receiver_coords,
466 lambda_mean,
467 mu_mean,
468 nthreads=nthreads,
469 rotate_sdn=int(rotate_sdn),
470 stack_sources=int(variant != 'normal'))
472 eps = 0.5 * (
473 results[:, :, 3:] +
474 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
476 dilatation \
477 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
479 stress_sdn = kron*lambda_mean*dilatation + 2.*mu_mean*eps
480 coefmat[:, idisl::3] = stress_sdn[:, :, (2, 5, 8)]\
481 .reshape(-1, npoints*3).T
482 else:
483 for isrc, source in enumerate(source_patches):
484 results = okada_ext.okada(
485 source.reshape(1, -1),
486 source_disl.reshape(1, -1),
487 receiver_coords,
488 lambda_mean,
489 mu_mean,
490 nthreads=nthreads,
491 rotate_sdn=int(rotate_sdn))
493 eps = 0.5 * (
494 results[:, 3:] +
495 results[:, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
497 dilatation \
498 = num.sum(eps[:, diag_ind], axis=1)[:, num.newaxis]
499 stress_sdn \
500 = kron * lambda_mean * dilatation+2. * mu_mean * eps
502 coefmat[:, isrc*3 + idisl] \
503 = stress_sdn[:, (2, 5, 8)].ravel()
505 if pure_shear:
506 coefmat[2::3, :] = 0.
508 return -coefmat / unit_disl
511def _make_okada_coefficient_matrix_slow(
512 source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1):
514 source_patches = num.array([
515 src.source_patch() for src in source_patches_list])
516 receiver_coords = source_patches[:, :3].copy()
518 npoints = len(source_patches_list)
520 if pure_shear:
521 n_eq = 2
522 else:
523 n_eq = 3
525 coefmat = num.zeros((npoints * 3, npoints * 3))
527 def ned2sdn_rotmat(strike, dip):
528 rotmat = mt.euler_to_matrix(
529 (dip + 180.) * d2r, strike * d2r, 0.)
530 return rotmat
532 lambda_mean = num.mean([src.lamb for src in source_patches_list])
533 shearmod_mean = num.mean([src.shearmod for src in source_patches_list])
535 unit_disl = 1.
536 disl_cases = {
537 'strikeslip': {
538 'slip': unit_disl,
539 'opening': 0.,
540 'rake': 0.},
541 'dipslip': {
542 'slip': unit_disl,
543 'opening': 0.,
544 'rake': 90.},
545 'tensileslip': {
546 'slip': 0.,
547 'opening': unit_disl,
548 'rake': 0.}
549 }
550 for idisl, case_type in enumerate([
551 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]):
552 case = disl_cases[case_type]
553 source_disl = num.array([
554 case['slip'] * num.cos(case['rake'] * d2r),
555 case['slip'] * num.sin(case['rake'] * d2r),
556 case['opening']])
558 for isource, source in enumerate(source_patches):
559 results = okada_ext.okada(
560 source[num.newaxis, :].copy(),
561 source_disl[num.newaxis, :].copy(),
562 receiver_coords,
563 lambda_mean,
564 shearmod_mean,
565 nthreads=nthreads,
566 rotate_sdn=int(rotate_sdn))
568 for irec in range(receiver_coords.shape[0]):
569 eps = num.zeros((3, 3))
570 for m in range(3):
571 for n in range(3):
572 eps[m, n] = 0.5 * (
573 results[irec][m * 3 + n + 3] +
574 results[irec][n * 3 + m + 3])
576 stress = num.zeros((3, 3))
577 dilatation = num.sum([eps[i, i] for i in range(3)])
579 for m, n in zip([0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]):
580 if m == n:
581 stress[m, n] = \
582 lambda_mean * \
583 dilatation + \
584 2. * shearmod_mean * \
585 eps[m, n]
587 else:
588 stress[m, n] = \
589 2. * shearmod_mean * \
590 eps[m, n]
591 stress[n, m] = stress[m, n]
593 normal = num.array([0., 0., -1.])
594 for isig in range(3):
595 tension = num.sum(stress[isig, :] * normal)
596 coefmat[irec * n_eq + isig, isource * n_eq + idisl] = \
597 tension / unit_disl
599 return coefmat
602def invert_fault_dislocations_bem(
603 stress_field,
604 coef_mat=None,
605 source_list=None,
606 pure_shear=False,
607 epsilon=None,
608 nthreads=1,
609 **kwargs):
610 '''
611 BEM least squares inversion to get fault dislocations given stress field.
613 Follows least squares inversion approach by Menke (1989) to calculate
614 dislocations on a fault with several segments from a given stress field.
615 The coefficient matrix connecting stresses and displacements of the fault
616 patches can either be specified by the user (``coef_mat``) or it is
617 calculated using the solution of Okada (1992) for a rectangular fault in a
618 homogeneous half space (``source_list``).
620 :param stress_field:
621 Stress change [Pa] for each source patch (as
622 ``stress_field[isource, icomponent]`` where isource indexes the source
623 patch and ``icomponent`` indexes component, ordered (strike, dip,
624 tensile).
625 :type stress_field:
626 :py:class:`~numpy.ndarray`: ``(nsources, 3)``
628 :param coef_mat:
629 Coefficient matrix connecting source patch dislocations and the stress
630 field.
631 :type coef_mat:
632 :py:class:`~numpy.ndarray`:
633 ``(len(source_list) * 3, len(source_list) * 3)``
635 :param source_list:
636 Source patches to be used for BEM.
637 :type source_list:
638 list of
639 :py:class:`~pyrocko.modelling.okada.OkadaSource`
641 :param epsilon:
642 If given, values in ``coef_mat`` smaller than ``epsilon`` are set to
643 zero.
644 :type epsilon:
645 float
647 :param nthreads:
648 Number of threads allowed.
649 :type nthreads:
650 int
652 :return:
653 Inverted displacements as ``displacements[isource, icomponent]``
654 where isource indexes the source patch and ``icomponent`` indexes
655 component, ordered (strike, dip, tensile).
656 :rtype:
657 :py:class:`~numpy.ndarray`: ``(nsources, 3)``
658 '''
660 if source_list is not None and coef_mat is None:
661 coef_mat = make_okada_coefficient_matrix(
662 source_list, pure_shear=pure_shear, nthreads=nthreads,
663 **kwargs)
665 if epsilon is not None:
666 coef_mat[coef_mat < epsilon] = 0.
668 idx = num.arange(0, coef_mat.shape[0])
669 if pure_shear:
670 idx = idx[idx % 3 != 2]
672 coef_mat_in = coef_mat[idx, :][:, idx]
673 disloc_est = num.zeros(coef_mat.shape[0])
675 if stress_field.ndim == 2:
676 stress_field = stress_field.ravel()
678 threadpool_limits = get_threadpool_limits()
680 with threadpool_limits(limits=nthreads, user_api='blas'):
681 try:
682 disloc_est[idx] = num.linalg.multi_dot([
683 num.linalg.inv(num.dot(coef_mat_in.T, coef_mat_in)),
684 coef_mat_in.T,
685 stress_field[idx]])
686 except num.linalg.LinAlgError as e:
687 logger.warning('Linear inversion failed!')
688 logger.warning(
689 'coef_mat: %s\nstress_field: %s',
690 coef_mat_in, stress_field[idx])
691 raise e
692 return disloc_est.reshape(-1, 3)
695__all__ = [
696 'AnalyticalSource',
697 'AnalyticalRectangularSource',
698 'OkadaSource',
699 'OkadaPatch',
700 'make_okada_coefficient_matrix',
701 'invert_fault_dislocations_bem']