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

5 

6''' 

7Elastostatic solutions and boundary element modelling for rectangular 

8dislocation sources. 

9''' 

10 

11import numpy as num 

12import logging 

13 

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 

19 

20guts_prefix = 'modelling' 

21 

22logger = logging.getLogger(__name__) 

23 

24d2r = num.pi/180. 

25r2d = 180./num.pi 

26km = 1e3 

27 

28 

29class AnalyticalSource(Location): 

30 ''' 

31 Base class for analytical source models. 

32 ''' 

33 

34 name = String.T( 

35 optional=True, 

36 default='') 

37 

38 time = Timestamp.T( 

39 default=0., 

40 help='Source origin time', 

41 optional=True) 

42 

43 vr = Float.T( 

44 default=0., 

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

46 optional=True) 

47 

48 @property 

49 def northing(self): 

50 return self.north_shift 

51 

52 @property 

53 def easting(self): 

54 return self.east_shift 

55 

56 

57class AnalyticalRectangularSource(AnalyticalSource): 

58 ''' 

59 Rectangular analytical source model. 

60 

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

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

63 ''' 

64 

65 strike = Float.T( 

66 default=0.0, 

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

68 

69 dip = Float.T( 

70 default=90.0, 

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

72 

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

77 

78 al1 = Float.T( 

79 default=0., 

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

81 

82 al2 = Float.T( 

83 default=0., 

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

85 

86 aw1 = Float.T( 

87 default=0., 

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

89 

90 aw2 = Float.T( 

91 default=0., 

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

93 

94 slip = Float.T( 

95 default=0., 

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

97 optional=True) 

98 

99 @property 

100 def length(self): 

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

102 

103 @property 

104 def width(self): 

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

106 

107 @property 

108 def area(self): 

109 return self.width * self.length 

110 

111 

112class OkadaSource(AnalyticalRectangularSource): 

113 ''' 

114 Rectangular Okada source model. 

115 ''' 

116 

117 opening = Float.T( 

118 default=0., 

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

120 optional=True) 

121 

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) 

130 

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) 

139 

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) 

149 

150 @property 

151 def poisson(self): 

152 if self.poisson__ is not None: 

153 return self.poisson__ 

154 

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

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

157 

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

159 

160 @poisson.setter 

161 def poisson(self, poisson): 

162 self.poisson__ = poisson 

163 

164 @property 

165 def lamb(self): 

166 

167 if self.lamb__ is not None: 

168 return self.lamb__ 

169 

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

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

172 

173 return ( 

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

175 

176 @lamb.setter 

177 def lamb(self, lamb): 

178 self.lamb__ = lamb 

179 

180 @property 

181 def shearmod(self): 

182 

183 if self.shearmod__ is not None: 

184 return self.shearmod__ 

185 

186 if self.poisson__ is None: 

187 raise ValueError('Poisson ratio is needed') 

188 

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

190 

191 @shearmod.setter 

192 def shearmod(self, shearmod): 

193 self.shearmod__ = shearmod 

194 

195 @property 

196 def seismic_moment(self): 

197 ''' 

198 Scalar Seismic moment :math:`M_0`. 

199 

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

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

202 

203 .. important :: 

204 

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

206 

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

209 

210 :return: 

211 Seismic moment release. 

212 :rtype: 

213 float 

214 ''' 

215 

216 mu = self.shearmod 

217 

218 disl = 0. 

219 if self.slip: 

220 disl = self.slip 

221 if self.opening: 

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

223 

224 return mu * self.area * disl 

225 

226 @property 

227 def moment_magnitude(self): 

228 ''' 

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

230 

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

232 

233 :returns: 

234 Moment magnitude. 

235 :rtype: 

236 float 

237 ''' 

238 return mt.moment_to_magnitude(self.seismic_moment) 

239 

240 def source_patch(self): 

241 ''' 

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

243 

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

245 

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

263 

264 def source_disloc(self): 

265 ''' 

266 Get source dislocation array for okada_ext.okada input. 

267 

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

269 source rake. 

270 

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

281 

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

283 ''' 

284 Discretize fault into rectilinear grid of fault patches. 

285 

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

287 sub-faults unchanged. 

288 

289 :param nlength: 

290 Number of patches in strike direction. 

291 :type nlength: 

292 int 

293 

294 :param nwidth: 

295 Number of patches in down-dip direction. 

296 :type nwidth: 

297 int 

298 

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 

306 

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

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

309 

310 patch_length = self.length / nlength 

311 patch_width = self.width / nwidth 

312 

313 al1 = -patch_length / 2. 

314 al2 = patch_length / 2. 

315 aw1 = -patch_width / 2. 

316 aw2 = patch_width / 2. 

317 

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. 

321 

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

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

324 

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

326 

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 

331 

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

337 

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) 

349 

350 

351class OkadaPatch(OkadaSource): 

352 

353 ''' 

354 Okada source with additional 2D indexes for bookkeeping. 

355 ''' 

356 

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

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

359 

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

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

362 self.parent = parent 

363 

364 

365def make_okada_coefficient_matrix( 

366 source_patches_list, 

367 pure_shear=False, 

368 rotate_sdn=True, 

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

370 

371 ''' 

372 Build coefficient matrix for given fault patches. 

373 

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. 

380 

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 

386 :param pure_shear: 

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

388 :type pure_shear: 

389 bool 

390 

391 :param rotate_sdn: 

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

393 :type rotate_sdn: 

394 bool 

395 

396 :param nthreads: 

397 Number of threads. 

398 :type nthreads: 

399 int 

400 

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

407 

408 if variant == 'slow': 

409 return _make_okada_coefficient_matrix_slow( 

410 source_patches_list, pure_shear, rotate_sdn, nthreads) 

411 

412 source_patches = num.array([ 

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

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

415 

416 npoints = len(source_patches_list) 

417 

418 if pure_shear: 

419 n_eq = 2 

420 else: 

421 n_eq = 3 

422 

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

424 

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

427 

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 } 

443 

444 diag_ind = [0, 4, 8] 

445 kron = num.zeros(9) 

446 kron[diag_ind] = 1. 

447 

448 if variant == 'normal': 

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

450 else: 

451 kron = kron[num.newaxis, :] 

452 

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

460 

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

471 

472 eps = 0.5 * ( 

473 results[:, :, 3:] + 

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

475 

476 dilatation \ 

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

478 

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

492 

493 eps = 0.5 * ( 

494 results[:, 3:] + 

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

496 

497 dilatation \ 

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

499 stress_sdn \ 

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

501 

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

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

504 

505 if pure_shear: 

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

507 

508 return -coefmat / unit_disl 

509 

510 

511def _make_okada_coefficient_matrix_slow( 

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

513 

514 source_patches = num.array([ 

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

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

517 

518 npoints = len(source_patches_list) 

519 

520 if pure_shear: 

521 n_eq = 2 

522 else: 

523 n_eq = 3 

524 

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

526 

527 def ned2sdn_rotmat(strike, dip): 

528 rotmat = mt.euler_to_matrix( 

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

530 return rotmat 

531 

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

534 

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

557 

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

567 

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

575 

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

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

578 

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] 

586 

587 else: 

588 stress[m, n] = \ 

589 2. * shearmod_mean * \ 

590 eps[m, n] 

591 stress[n, m] = stress[m, n] 

592 

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 

598 

599 return coefmat 

600 

601 

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. 

612 

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

619 

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

627 

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

634 

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` 

640 

641 :param epsilon: 

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

643 zero. 

644 :type epsilon: 

645 float 

646 

647 :param nthreads: 

648 Number of threads allowed. 

649 :type nthreads: 

650 int 

651 

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

659 

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) 

664 

665 if epsilon is not None: 

666 coef_mat[coef_mat < epsilon] = 0. 

667 

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

669 if pure_shear: 

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

671 

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

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

674 

675 if stress_field.ndim == 2: 

676 stress_field = stress_field.ravel() 

677 

678 threadpool_limits = get_threadpool_limits() 

679 

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) 

693 

694 

695__all__ = [ 

696 'AnalyticalSource', 

697 'AnalyticalRectangularSource', 

698 'OkadaSource', 

699 'OkadaPatch', 

700 'make_okada_coefficient_matrix', 

701 'invert_fault_dislocations_bem']