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

125 

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) 

134 

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) 

144 

145 @property 

146 def poisson(self): 

147 if self.poisson__ is not None: 

148 return self.poisson__ 

149 

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

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

152 

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

154 

155 @poisson.setter 

156 def poisson(self, poisson): 

157 self.poisson__ = poisson 

158 

159 @property 

160 def lamb(self): 

161 

162 if self.lamb__ is not None: 

163 return self.lamb__ 

164 

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

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

167 

168 return ( 

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

170 

171 @lamb.setter 

172 def lamb(self, lamb): 

173 self.lamb__ = lamb 

174 

175 @property 

176 def shearmod(self): 

177 

178 if self.shearmod__ is not None: 

179 return self.shearmod__ 

180 

181 if self.poisson__ is None: 

182 raise ValueError('Poisson ratio is needed') 

183 

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

185 

186 @shearmod.setter 

187 def shearmod(self, shearmod): 

188 self.shearmod__ = shearmod 

189 

190 @property 

191 def seismic_moment(self): 

192 ''' 

193 Scalar Seismic moment :math:`M_0`. 

194 

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

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

197 

198 .. important :: 

199 

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

201 

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

204 

205 :return: 

206 Seismic moment release. 

207 :rtype: 

208 float 

209 ''' 

210 

211 mu = self.shearmod 

212 

213 disl = 0. 

214 if self.slip: 

215 disl = self.slip 

216 if self.opening: 

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

218 

219 return mu * self.area * disl 

220 

221 @property 

222 def moment_magnitude(self): 

223 ''' 

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

225 

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

227 

228 :returns: 

229 Moment magnitude. 

230 :rtype: 

231 float 

232 ''' 

233 return mt.moment_to_magnitude(self.seismic_moment) 

234 

235 def source_patch(self): 

236 ''' 

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

238 

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

240 

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

258 

259 def source_disloc(self): 

260 ''' 

261 Get source dislocation array for okada_ext.okada input. 

262 

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

264 source rake. 

265 

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

276 

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

278 ''' 

279 Discretize fault into rectilinear grid of fault patches. 

280 

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

282 sub-faults unchanged. 

283 

284 :param nlength: 

285 Number of patches in strike direction. 

286 :type nlength: 

287 int 

288 

289 :param nwidth: 

290 Number of patches in down-dip direction. 

291 :type nwidth: 

292 int 

293 

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 

301 

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

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

304 

305 patch_length = self.length / nlength 

306 patch_width = self.width / nwidth 

307 

308 al1 = -patch_length / 2. 

309 al2 = patch_length / 2. 

310 aw1 = -patch_width / 2. 

311 aw2 = patch_width / 2. 

312 

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. 

316 

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

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

319 

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

321 

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 

326 

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

332 

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) 

344 

345 

346class OkadaPatch(OkadaSource): 

347 

348 ''' 

349 Okada source with additional 2D indexes for bookkeeping. 

350 ''' 

351 

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

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

354 

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

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

357 self.parent = parent 

358 

359 

360def make_okada_coefficient_matrix( 

361 source_patches_list, 

362 pure_shear=False, 

363 rotate_sdn=True, 

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

365 

366 ''' 

367 Build coefficient matrix for given fault patches. 

368 

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. 

375 

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

380 

381 :param pure_shear: 

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

383 :type pure_shear: 

384 optional, bool 

385 

386 :param rotate_sdn: 

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

388 :type rotate_sdn: 

389 optional, bool 

390 

391 :param nthreads: 

392 Number of threads. 

393 :type nthreads: 

394 optional, int 

395 

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

402 

403 if variant == 'slow': 

404 return _make_okada_coefficient_matrix_slow( 

405 source_patches_list, pure_shear, rotate_sdn, nthreads) 

406 

407 source_patches = num.array([ 

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

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

410 

411 npoints = len(source_patches_list) 

412 

413 if pure_shear: 

414 n_eq = 2 

415 else: 

416 n_eq = 3 

417 

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

419 

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

422 

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 } 

438 

439 diag_ind = [0, 4, 8] 

440 kron = num.zeros(9) 

441 kron[diag_ind] = 1. 

442 

443 if variant == 'normal': 

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

445 else: 

446 kron = kron[num.newaxis, :] 

447 

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

455 

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

466 

467 eps = 0.5 * ( 

468 results[:, :, 3:] + 

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

470 

471 dilatation \ 

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

473 

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

487 

488 eps = 0.5 * ( 

489 results[:, 3:] + 

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

491 

492 dilatation \ 

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

494 stress_sdn \ 

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

496 

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

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

499 

500 if pure_shear: 

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

502 

503 return -coefmat / unit_disl 

504 

505 

506def _make_okada_coefficient_matrix_slow( 

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

508 

509 source_patches = num.array([ 

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

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

512 

513 npoints = len(source_patches_list) 

514 

515 if pure_shear: 

516 n_eq = 2 

517 else: 

518 n_eq = 3 

519 

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

521 

522 def ned2sdn_rotmat(strike, dip): 

523 rotmat = mt.euler_to_matrix( 

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

525 return rotmat 

526 

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

529 

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

552 

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

562 

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

570 

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

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

573 

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] 

581 

582 else: 

583 stress[m, n] = \ 

584 2. * shearmod_mean * \ 

585 eps[m, n] 

586 stress[n, m] = stress[m, n] 

587 

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 

593 

594 return coefmat 

595 

596 

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. 

607 

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

614 

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

622 

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

629 

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` 

635 

636 :param epsilon: 

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

638 zero. 

639 :type epsilon: 

640 optional, float 

641 

642 :param nthreads: 

643 Number of threads allowed. 

644 :type nthreads: 

645 int 

646 

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

654 

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) 

659 

660 if epsilon is not None: 

661 coef_mat[coef_mat < epsilon] = 0. 

662 

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

664 if pure_shear: 

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

666 

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

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

669 

670 if stress_field.ndim == 2: 

671 stress_field = stress_field.ravel() 

672 

673 threadpool_limits = get_threadpool_limits() 

674 

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) 

688 

689 

690__all__ = [ 

691 'AnalyticalSource', 

692 'AnalyticalRectangularSource', 

693 'OkadaSource', 

694 'OkadaPatch', 

695 'make_okada_coefficient_matrix', 

696 'invert_fault_dislocations_bem']