1# https://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6import numpy as num 

7import logging 

8 

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 

14 

15guts_prefix = 'modelling' 

16 

17logger = logging.getLogger(__name__) 

18 

19d2r = num.pi/180. 

20r2d = 180./num.pi 

21km = 1e3 

22 

23 

24class AnalyticalSource(Location): 

25 ''' 

26 Base class for analytical source models. 

27 ''' 

28 

29 name = String.T( 

30 optional=True, 

31 default='') 

32 

33 time = Timestamp.T( 

34 default=0., 

35 help='Source origin time', 

36 optional=True) 

37 

38 vr = Float.T( 

39 default=0., 

40 help='Rupture velocity [m/s]', 

41 optional=True) 

42 

43 @property 

44 def northing(self): 

45 return self.north_shift 

46 

47 @property 

48 def easting(self): 

49 return self.east_shift 

50 

51 

52class AnalyticalRectangularSource(AnalyticalSource): 

53 ''' 

54 Rectangular analytical source model. 

55 

56 Coordinates on the source plane are with respect to the origin point given 

57 by `(lat, lon, east_shift, north_shift, depth)`. 

58 ''' 

59 

60 strike = Float.T( 

61 default=0.0, 

62 help='Strike direction in [deg], measured clockwise from north.') 

63 

64 dip = Float.T( 

65 default=90.0, 

66 help='Dip angle in [deg], measured downward from horizontal.') 

67 

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.') 

72 

73 al1 = Float.T( 

74 default=0., 

75 help='Left edge source plane coordinate [m].') 

76 

77 al2 = Float.T( 

78 default=0., 

79 help='Right edge source plane coordinate [m].') 

80 

81 aw1 = Float.T( 

82 default=0., 

83 help='Lower edge source plane coordinate [m].') 

84 

85 aw2 = Float.T( 

86 default=0., 

87 help='Upper edge source plane coordinate [m].') 

88 

89 slip = Float.T( 

90 default=0., 

91 help='Slip on the rectangular source area [m].', 

92 optional=True) 

93 

94 @property 

95 def length(self): 

96 return abs(-self.al1 + self.al2) 

97 

98 @property 

99 def width(self): 

100 return abs(-self.aw1 + self.aw2) 

101 

102 @property 

103 def area(self): 

104 return self.width * self.length 

105 

106 

107class OkadaSource(AnalyticalRectangularSource): 

108 ''' 

109 Rectangular Okada source model. 

110 ''' 

111 

112 opening = Float.T( 

113 default=0., 

114 help='Opening of the plane in [m].', 

115 optional=True) 

116 

117 poisson__ = Float.T( 

118 default=0.25, 

119 help='Poisson\'s ratio :math:`\\nu`.', 

120 optional=True) 

121 

122 lamb__ = Float.T( 

123 help='First Lame\' s parameter :math:`\\lambda` [Pa].', 

124 optional=True) 

125 

126 shearmod__ = Float.T( 

127 default=32.0e9, 

128 help='Shear modulus along the plane :math:`\\mu` [Pa].', 

129 optional=True) 

130 

131 @property 

132 def poisson(self): 

133 ''' 

134 Poisson\' s ratio :math:`\\nu` (if not given). 

135 

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 ''' 

140 

141 if self.poisson__ is not None: 

142 return self.poisson__ 

143 

144 if self.shearmod__ is None or self.lamb__ is None: 

145 raise ValueError('Shearmod and lambda are needed') 

146 

147 return (self.lamb__) / (2. * (self.lamb__ + self.shearmod__)) 

148 

149 @poisson.setter 

150 def poisson(self, poisson): 

151 self.poisson__ = poisson 

152 

153 @property 

154 def lamb(self): 

155 ''' 

156 First Lame\' s parameter :math:`\\lambda` (if not given). 

157 

158 Poisson\' s ratio :math:`\\nu` and shear modulus :math:`\\mu` must be 

159 available to calculate the first Lame\' s parameter :math:`\\lambda`. 

160 

161 .. important :: 

162 

163 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`. 

164 

165 Through :math:`\\nu = \\frac{\\lambda}{2(\\lambda + \\mu)}` this 

166 leads to :math:`\\lambda = \\frac{2 \\mu \\nu}{1-2\\nu}`. 

167 

168 ''' 

169 

170 if self.lamb__ is not None: 

171 return self.lamb__ 

172 

173 if self.shearmod__ is None or self.poisson__ is None: 

174 raise ValueError('Shearmod and poisson ratio are needed') 

175 

176 return ( 

177 2. * self.poisson__ * self.shearmod__) / (1. - 2. * self.poisson__) 

178 

179 @lamb.setter 

180 def lamb(self, lamb): 

181 self.lamb__ = lamb 

182 

183 @property 

184 def shearmod(self): 

185 ''' 

186 Shear modulus :math:`\\mu` (if not given). 

187 

188 Poisson ratio\' s :math:`\\nu` must be available. 

189 

190 .. important :: 

191 

192 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`. 

193 

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}`. 

196 

197 ''' 

198 

199 if self.shearmod__ is not None: 

200 return self.shearmod__ 

201 

202 if self.poisson__ is None: 

203 raise ValueError('Poisson ratio is needed') 

204 

205 return (8. * (1. + self.poisson__)) / (1. - 2. * self.poisson__) 

206 

207 @shearmod.setter 

208 def shearmod(self, shearmod): 

209 self.shearmod__ = shearmod 

210 

211 @property 

212 def seismic_moment(self): 

213 ''' 

214 Scalar Seismic moment :math:`M_0`. 

215 

216 Code copied from Kite. It disregards the opening (as for now). 

217 We assume :math:`M_0 = mu A D`. 

218 

219 .. important :: 

220 

221 We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`. 

222 

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}`. 

225 

226 :return: 

227 Seismic moment release. 

228 :rtype: 

229 float 

230 ''' 

231 

232 mu = self.shearmod 

233 

234 disl = 0. 

235 if self.slip: 

236 disl = self.slip 

237 if self.opening: 

238 disl = (disl**2 + self.opening**2)**.5 

239 

240 return mu * self.area * disl 

241 

242 @property 

243 def moment_magnitude(self): 

244 ''' 

245 Moment magnitude :math:`M_\\mathrm{w}` from seismic moment. 

246 

247 We assume :math:`M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0) - 10.7`. 

248 

249 :returns: 

250 Moment magnitude. 

251 :rtype: 

252 float 

253 ''' 

254 return mt.moment_to_magnitude(self.seismic_moment) 

255 

256 def source_patch(self): 

257 ''' 

258 Get source location and geometry array for okada_ext.okada input. 

259 

260 The values are defined according to Okada (1992). 

261 

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]) 

279 

280 def source_disloc(self): 

281 ''' 

282 Get source dislocation array for okada_ext.okada input. 

283 

284 The given slip is splitted into a strike and an updip part based on the 

285 source rake. 

286 

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]) 

297 

298 def discretize(self, nlength, nwidth, *args, **kwargs): 

299 ''' 

300 Discretize fault into rectilinear grid of fault patches. 

301 

302 Fault orientation, slip and elastic parameters are passed to the 

303 sub-faults unchanged. 

304 

305 :param nlength: 

306 Number of patches in strike direction. 

307 :type nlength: 

308 int 

309 

310 :param nwidth: 

311 Number of patches in down-dip direction. 

312 :type nwidth: 

313 int 

314 

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 

322 

323 il = num.repeat(num.arange(nlength), nwidth) 

324 iw = num.tile(num.arange(nwidth), nlength) 

325 

326 patch_length = self.length / nlength 

327 patch_width = self.width / nwidth 

328 

329 al1 = -patch_length / 2. 

330 al2 = patch_length / 2. 

331 aw1 = -patch_width / 2. 

332 aw2 = patch_width / 2. 

333 

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. 

337 

338 source_points[:, 0] += self.al1 

339 source_points[:, 1] -= self.aw2 

340 

341 rotmat = mt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.) 

342 

343 source_points_rot = num.dot(rotmat.T, source_points.T).T 

344 source_points_rot[:, 0] += self.northing 

345 source_points_rot[:, 1] += self.easting 

346 source_points_rot[:, 2] += self.depth 

347 

348 kwargs = { 

349 prop: getattr(self, prop) for prop in self.T.propnames 

350 if prop not in [ 

351 'north_shift', 'east_shift', 'depth', 

352 'al1', 'al2', 'aw1', 'aw2']} 

353 

354 return ( 

355 [OkadaPatch( 

356 parent=self, 

357 ix=src_point[0], 

358 iy=src_point[1], 

359 north_shift=coord[0], 

360 east_shift=coord[1], 

361 depth=coord[2], 

362 al1=al1, al2=al2, aw1=aw1, aw2=aw2, **kwargs) 

363 for src_point, coord in zip(source_points, source_points_rot)], 

364 source_points) 

365 

366 

367class OkadaPatch(OkadaSource): 

368 

369 ''' 

370 Okada source with additional 2D indexes for bookkeeping. 

371 ''' 

372 

373 ix = Int.T(help='Relative index of the patch in x') 

374 iy = Int.T(help='Relative index of the patch in y') 

375 

376 def __init__(self, parent=None, *args, **kwargs): 

377 OkadaSource.__init__(self, *args, **kwargs) 

378 self.parent = parent 

379 

380 

381def make_okada_coefficient_matrix( 

382 source_patches_list, 

383 pure_shear=False, 

384 rotate_sdn=True, 

385 nthreads=1, variant='normal'): 

386 

387 ''' 

388 Build coefficient matrix for given fault patches. 

389 

390 The boundary element method (BEM) for a discretized fault and the 

391 determination of the slip distribution :math:`\\Delta u` from stress drop 

392 :math:`\\Delta \\sigma` is based on 

393 :math:`\\Delta \\sigma = \\mathbf{C} \\cdot \\Delta u`. Here the 

394 coefficient matrix :math:`\\mathbf{C}` is built, based on the displacements 

395 from Okada's solution (Okada, 1992) and their partial derivatives. 

396 

397 :param source_patches_list: 

398 Source patches, to be used in BEM. 

399 :type source_patches_list: 

400 list of :py:class:`~pyrocko.modelling.okada.OkadaSource`. 

401 

402 :param pure_shear: 

403 If ``True``, only shear forces are taken into account. 

404 :type pure_shear: 

405 optional, bool 

406 

407 :param rotate_sdn: 

408 If ``True``, rotate to strike, dip, normal. 

409 :type rotate_sdn: 

410 optional, bool 

411 

412 :param nthreads: 

413 Number of threads. 

414 :type nthreads: 

415 optional, int 

416 

417 :return: 

418 Coefficient matrix for all source combinations. 

419 :rtype: 

420 :py:class:`~numpy.ndarray`: 

421 ``(len(source_patches_list) * 3, len(source_patches_list) * 3)`` 

422 ''' 

423 

424 if variant == 'slow': 

425 return _make_okada_coefficient_matrix_slow( 

426 source_patches_list, pure_shear, rotate_sdn, nthreads) 

427 

428 source_patches = num.array([ 

429 src.source_patch() for src in source_patches_list]) 

430 receiver_coords = source_patches[:, :3].copy() 

431 

432 npoints = len(source_patches_list) 

433 

434 if pure_shear: 

435 n_eq = 2 

436 else: 

437 n_eq = 3 

438 

439 coefmat = num.zeros((npoints * 3, npoints * 3)) 

440 

441 lambda_mean = num.mean([src.lamb for src in source_patches_list]) 

442 mu_mean = num.mean([src.shearmod for src in source_patches_list]) 

443 

444 unit_disl = 1. 

445 disl_cases = { 

446 'strikeslip': { 

447 'slip': unit_disl, 

448 'opening': 0., 

449 'rake': 0.}, 

450 'dipslip': { 

451 'slip': unit_disl, 

452 'opening': 0., 

453 'rake': 90.}, 

454 'tensileslip': { 

455 'slip': 0., 

456 'opening': unit_disl, 

457 'rake': 0.} 

458 } 

459 

460 diag_ind = [0, 4, 8] 

461 kron = num.zeros(9) 

462 kron[diag_ind] = 1. 

463 

464 if variant == 'normal': 

465 kron = kron[num.newaxis, num.newaxis, :] 

466 else: 

467 kron = kron[num.newaxis, :] 

468 

469 for idisl, case_type in enumerate([ 

470 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]): 

471 case = disl_cases[case_type] 

472 source_disl = num.array([ 

473 case['slip'] * num.cos(case['rake'] * d2r), 

474 case['slip'] * num.sin(case['rake'] * d2r), 

475 case['opening']]) 

476 

477 if variant == 'normal': 

478 results = okada_ext.okada( 

479 source_patches, 

480 num.tile(source_disl, npoints).reshape(-1, 3), 

481 receiver_coords, 

482 lambda_mean, 

483 mu_mean, 

484 nthreads=nthreads, 

485 rotate_sdn=int(rotate_sdn), 

486 stack_sources=int(variant != 'normal')) 

487 

488 eps = 0.5 * ( 

489 results[:, :, 3:] + 

490 results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)]) 

491 

492 dilatation \ 

493 = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis] 

494 

495 stress_sdn = kron*lambda_mean*dilatation + 2.*mu_mean*eps 

496 coefmat[:, idisl::3] = stress_sdn[:, :, (2, 5, 8)]\ 

497 .reshape(-1, npoints*3).T 

498 else: 

499 for isrc, source in enumerate(source_patches): 

500 results = okada_ext.okada( 

501 source[num.newaxis, :], 

502 source_disl[num.newaxis, :], 

503 receiver_coords, 

504 lambda_mean, 

505 mu_mean, 

506 nthreads=nthreads, 

507 rotate_sdn=int(rotate_sdn)) 

508 

509 eps = 0.5 * ( 

510 results[:, 3:] + 

511 results[:, (3, 6, 9, 4, 7, 10, 5, 8, 11)]) 

512 

513 dilatation \ 

514 = num.sum(eps[:, diag_ind], axis=1)[:, num.newaxis] 

515 stress_sdn \ 

516 = kron * lambda_mean * dilatation+2. * mu_mean * eps 

517 

518 coefmat[:, isrc*3 + idisl] \ 

519 = stress_sdn[:, (2, 5, 8)].ravel() 

520 

521 if pure_shear: 

522 coefmat[2::3, :] = 0. 

523 

524 return -coefmat / unit_disl 

525 

526 

527def _make_okada_coefficient_matrix_slow( 

528 source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1): 

529 

530 source_patches = num.array([ 

531 src.source_patch() for src in source_patches_list]) 

532 receiver_coords = source_patches[:, :3].copy() 

533 

534 npoints = len(source_patches_list) 

535 

536 if pure_shear: 

537 n_eq = 2 

538 else: 

539 n_eq = 3 

540 

541 coefmat = num.zeros((npoints * 3, npoints * 3)) 

542 

543 def ned2sdn_rotmat(strike, dip): 

544 rotmat = mt.euler_to_matrix( 

545 (dip + 180.) * d2r, strike * d2r, 0.) 

546 return rotmat 

547 

548 lambda_mean = num.mean([src.lamb for src in source_patches_list]) 

549 shearmod_mean = num.mean([src.shearmod for src in source_patches_list]) 

550 

551 unit_disl = 1. 

552 disl_cases = { 

553 'strikeslip': { 

554 'slip': unit_disl, 

555 'opening': 0., 

556 'rake': 0.}, 

557 'dipslip': { 

558 'slip': unit_disl, 

559 'opening': 0., 

560 'rake': 90.}, 

561 'tensileslip': { 

562 'slip': 0., 

563 'opening': unit_disl, 

564 'rake': 0.} 

565 } 

566 for idisl, case_type in enumerate([ 

567 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]): 

568 case = disl_cases[case_type] 

569 source_disl = num.array([ 

570 case['slip'] * num.cos(case['rake'] * d2r), 

571 case['slip'] * num.sin(case['rake'] * d2r), 

572 case['opening']]) 

573 

574 for isource, source in enumerate(source_patches): 

575 results = okada_ext.okada( 

576 source[num.newaxis, :].copy(), 

577 source_disl[num.newaxis, :].copy(), 

578 receiver_coords, 

579 lambda_mean, 

580 shearmod_mean, 

581 nthreads=nthreads, 

582 rotate_sdn=int(rotate_sdn)) 

583 

584 for irec in range(receiver_coords.shape[0]): 

585 eps = num.zeros((3, 3)) 

586 for m in range(3): 

587 for n in range(3): 

588 eps[m, n] = 0.5 * ( 

589 results[irec][m * 3 + n + 3] + 

590 results[irec][n * 3 + m + 3]) 

591 

592 stress = num.zeros((3, 3)) 

593 dilatation = num.sum([eps[i, i] for i in range(3)]) 

594 

595 for m, n in zip([0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]): 

596 if m == n: 

597 stress[m, n] = \ 

598 lambda_mean * \ 

599 dilatation + \ 

600 2. * shearmod_mean * \ 

601 eps[m, n] 

602 

603 else: 

604 stress[m, n] = \ 

605 2. * shearmod_mean * \ 

606 eps[m, n] 

607 stress[n, m] = stress[m, n] 

608 

609 normal = num.array([0., 0., -1.]) 

610 for isig in range(3): 

611 tension = num.sum(stress[isig, :] * normal) 

612 coefmat[irec * n_eq + isig, isource * n_eq + idisl] = \ 

613 tension / unit_disl 

614 

615 return coefmat 

616 

617 

618def invert_fault_dislocations_bem( 

619 stress_field, 

620 coef_mat=None, 

621 source_list=None, 

622 pure_shear=False, 

623 epsilon=None, 

624 nthreads=1, 

625 **kwargs): 

626 ''' 

627 BEM least squares inversion to get fault dislocations given stress field. 

628 

629 Follows least squares inversion approach by Menke (1989) to calculate 

630 dislocations on a fault with several segments from a given stress field. 

631 The coefficient matrix connecting stresses and displacements of the fault 

632 patches can either be specified by the user (``coef_mat``) or it is 

633 calculated using the solution of Okada (1992) for a rectangular fault in a 

634 homogeneous half space (``source_list``). 

635 

636 :param stress_field: 

637 Stress change [Pa] for each source patch (as 

638 ``stress_field[isource, icomponent]`` where isource indexes the source 

639 patch and ``icomponent`` indexes component, ordered (strike, dip, 

640 tensile). 

641 :type stress_field: 

642 :py:class:`~numpy.ndarray`: ``(nsources, 3)`` 

643 

644 :param coef_mat: 

645 Coefficient matrix connecting source patch dislocations and the stress 

646 field. 

647 :type coef_mat: 

648 optional, :py:class:`~numpy.ndarray`: 

649 ``(len(source_list) * 3, len(source_list) * 3)`` 

650 

651 :param source_list: 

652 Source patches to be used for BEM. 

653 :type source_list: 

654 optional, list of 

655 :py:class:`~pyrocko.modelling.okada.OkadaSource` 

656 

657 :param epsilon: 

658 If given, values in ``coef_mat`` smaller than ``epsilon`` are set to 

659 zero. 

660 :type epsilon: 

661 optional, float 

662 

663 :param nthreads: 

664 Number of threads allowed. 

665 :type nthreads: 

666 int 

667 

668 :return: 

669 Inverted displacements as ``displacements[isource, icomponent]`` 

670 where isource indexes the source patch and ``icomponent`` indexes 

671 component, ordered (strike, dip, tensile). 

672 :rtype: 

673 :py:class:`~numpy.ndarray`: ``(nsources, 3)`` 

674 ''' 

675 

676 if source_list is not None and coef_mat is None: 

677 coef_mat = make_okada_coefficient_matrix( 

678 source_list, pure_shear=pure_shear, nthreads=nthreads, 

679 **kwargs) 

680 

681 if epsilon is not None: 

682 coef_mat[coef_mat < epsilon] = 0. 

683 

684 idx = num.arange(0, coef_mat.shape[0]) 

685 if pure_shear: 

686 idx = idx[idx % 3 != 2] 

687 

688 coef_mat_in = coef_mat[idx, :][:, idx] 

689 disloc_est = num.zeros(coef_mat.shape[0]) 

690 

691 if stress_field.ndim == 2: 

692 stress_field = stress_field.ravel() 

693 

694 threadpool_limits = get_threadpool_limits() 

695 

696 with threadpool_limits(limits=nthreads, user_api='blas'): 

697 try: 

698 disloc_est[idx] = num.linalg.multi_dot([ 

699 num.linalg.inv(num.dot(coef_mat_in.T, coef_mat_in)), 

700 coef_mat_in.T, 

701 stress_field[idx]]) 

702 except num.linalg.LinAlgError as e: 

703 logger.warning('Linear inversion failed!') 

704 logger.warning( 

705 'coef_mat: %s\nstress_field: %s', 

706 coef_mat_in, stress_field[idx]) 

707 raise e 

708 return disloc_est.reshape(-1, 3) 

709 

710 

711__all__ = [ 

712 'AnalyticalSource', 

713 'AnalyticalRectangularSource', 

714 'OkadaSource', 

715 'OkadaPatch', 

716 'make_okada_coefficient_matrix', 

717 'invert_fault_dislocations_bem']