Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gf/tractions.py: 70%

166 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 06:59 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Traction field parameterizations used by :doc:`PDR source 

8</topics/pseudo-dynamic-rupture>`. 

9''' 

10 

11import logging 

12import numpy as num 

13from pyrocko.guts import Object, Float, List, StringChoice, Int 

14from pyrocko.guts_array import Array 

15 

16logger = logging.getLogger('pyrocko.gf.tractions') 

17km = 1e3 

18d2r = num.pi/180. 

19r2d = 180./num.pi 

20 

21 

22def tukey_window(N, alpha): 

23 assert alpha <= 1. 

24 window = num.ones(N) 

25 n = num.arange(N) 

26 

27 N_f = int((alpha * N)//2) 

28 window[:N_f] = .5 * (1. - num.cos((2*num.pi * n[:N_f])/(alpha * N))) 

29 window[(N-N_f):] = window[:N_f][::-1] 

30 return window 

31 

32 

33def planck_window(N, epsilon): 

34 assert epsilon <= 1. 

35 window = num.ones(N) 

36 n = num.arange(N) 

37 

38 N_f = int((epsilon * N)) 

39 window[:N_f] = \ 

40 (1. + num.exp((epsilon * N) / n[:N_f] - 

41 (epsilon * N) / ((epsilon * N - n[:N_f]))))**-1. 

42 window[(N-N_f):] = window[:N_f][::-1] 

43 return window 

44 

45 

46class AbstractTractionField(Object): 

47 ''' 

48 Base class for multiplicative traction fields (tapers). 

49 

50 Fields of this type a re multiplied in the 

51 :py:class:`~pyrocko.gf.tractions.TractionComposition` 

52 ''' 

53 operation = 'mult' 

54 

55 def get_tractions(self, nx, ny, patches): 

56 raise NotImplementedError 

57 

58 

59class TractionField(AbstractTractionField): 

60 ''' 

61 Base class for additive traction fields. 

62 

63 Fields of this type are added in the 

64 :py:class:`~pyrocko.gf.tractions.TractionComposition` 

65 ''' 

66 operation = 'add' 

67 

68 def get_tractions(self, nx, ny, patches): 

69 raise NotImplementedError 

70 

71 

72class TractionComposition(TractionField): 

73 ''' 

74 Composition of traction fields. 

75 

76 :py:class:`~pyrocko.gf.tractions.TractionField` and 

77 :py:class:`~pyrocko.gf.tractions.AbstractTractionField` can be combined 

78 to realize a combination of different fields. 

79 ''' 

80 components = List.T( 

81 AbstractTractionField.T(), 

82 default=[], 

83 help='Ordered list of tractions.') 

84 

85 def get_tractions(self, nx, ny, patches=None): 

86 npatches = nx * ny 

87 tractions = num.zeros((npatches, 3)) 

88 

89 for comp in self.components: 

90 if comp.operation == 'add': 

91 tractions += comp.get_tractions(nx, ny, patches) 

92 elif comp.operation == 'mult': 

93 tractions *= comp.get_tractions(nx, ny, patches) 

94 else: 

95 raise AttributeError( 

96 'Component %s has an invalid operation %s.' % 

97 (comp, comp.operation)) 

98 

99 return tractions 

100 

101 def add_component(self, field): 

102 logger.debug('Adding traction component.') 

103 self.components.append(field) 

104 

105 

106class HomogeneousTractions(TractionField): 

107 ''' 

108 Homogeneous traction field. 

109 

110 The traction vectors in strike, dip and normal direction are acting 

111 homogeneously on the rupture plane. 

112 ''' 

113 

114 strike = Float.T( 

115 default=1., 

116 help='Tractions in strike direction [Pa].') 

117 dip = Float.T( 

118 default=1., 

119 help='Traction in dip direction (up) [Pa].') 

120 normal = Float.T( 

121 default=1., 

122 help='Traction in normal direction [Pa].') 

123 

124 def get_tractions(self, nx, ny, patches=None): 

125 npatches = nx * ny 

126 

127 return num.tile( 

128 (self.strike, self.dip, self.normal), npatches) \ 

129 .reshape(-1, 3) 

130 

131 

132class DirectedTractions(TractionField): 

133 ''' 

134 Directed traction field. 

135 

136 The traction vectors are following a uniform ``rake``. 

137 ''' 

138 

139 rake = Float.T( 

140 default=0., 

141 help='Rake angle in [deg], ' 

142 'measured counter-clockwise from right-horizontal ' 

143 'in on-plane view. Rake is translated into homogenous tractions ' 

144 'in strike and up-dip direction.') 

145 traction = Float.T( 

146 default=1., 

147 help='Traction in rake direction [Pa].') 

148 

149 def get_tractions(self, nx, ny, patches=None): 

150 npatches = nx * ny 

151 

152 strike = num.cos(self.rake*d2r) * self.traction 

153 dip = num.sin(self.rake*d2r) * self.traction 

154 normal = 0. 

155 

156 return num.tile((strike, dip, normal), npatches).reshape(-1, 3) 

157 

158 

159class SelfSimilarTractions(TractionField): 

160 ''' 

161 Traction model following Power & Tullis (1991). 

162 

163 The traction vectors are calculated as a sum of 2D-cosines with a constant 

164 amplitude / wavelength ratio. The wavenumber kx and ky are constant for 

165 each cosine function. The rank defines the maximum wavenumber used for 

166 summation. So, e.g. a rank of 3 will lead to a summation of cosines with 

167 ``kx = ky`` in (1, 2, 3). 

168 Each cosine has an associated phases, which defines both the phase shift 

169 and also the shift from the rupture plane centre. 

170 Finally the summed cosines are translated into shear tractions based on the 

171 rake and normalized with ``traction_max``. 

172 

173 ''' 

174 rank = Int.T( 

175 default=1, 

176 help='Maximum summed cosine wavenumber/spatial frequency.') 

177 

178 rake = Float.T( 

179 default=0., 

180 help='Rake angle in [deg], ' 

181 'measured counter-clockwise from right-horizontal ' 

182 'in on-plane view. Rake is translated into homogenous tractions ' 

183 'in strike and up-dip direction.') 

184 

185 traction_max = Float.T( 

186 default=1., 

187 help='Maximum traction vector length [Pa].') 

188 

189 phases = Array.T( 

190 optional=True, 

191 dtype=num.float64, 

192 shape=(None,), 

193 help='Phase shift of the cosines in [rad].') 

194 

195 def get_phases(self): 

196 if self.phases is not None: 

197 if self.phases.shape[0] == self.rank: 

198 return self.phases 

199 

200 return (num.random.random(self.rank) * 2. - 1.) * num.pi 

201 

202 def get_tractions(self, nx, ny, patches=None): 

203 z = num.zeros((ny, nx)) 

204 phases = self.get_phases() 

205 

206 for i in range(1, self.rank+1): 

207 x = num.linspace(-i*num.pi, i*num.pi, nx) + i*phases[i-1] 

208 y = num.linspace(-i*num.pi, i*num.pi, ny) + i*phases[i-1] 

209 x, y = num.meshgrid(x, y) 

210 r = num.sqrt(x**2 + y**2) 

211 z += 1. / i * num.cos(r + phases[i-1]) 

212 

213 t = num.zeros((nx*ny, 3)) 

214 t[:, 0] = num.cos(self.rake*d2r) * z.ravel(order='F') 

215 t[:, 1] = num.sin(self.rake*d2r) * z.ravel(order='F') 

216 

217 t *= self.traction_max / num.max(num.linalg.norm(t, axis=1)) 

218 

219 return t 

220 

221 

222class FractalTractions(TractionField): 

223 ''' 

224 Fractal traction field. 

225 ''' 

226 

227 rseed = Int.T( 

228 default=None, 

229 optional=True, 

230 help='Seed for :py:class:`~numpy.random.RandomState`.' 

231 'If ``None``, an random seed will be initialized.') 

232 

233 rake = Float.T( 

234 default=0., 

235 help='Rake angle in [deg], ' 

236 'measured counter-clockwise from right-horizontal ' 

237 'in on-plane view. Rake is translated into homogenous tractions ' 

238 'in strike and up-dip direction.') 

239 

240 traction_max = Float.T( 

241 default=1., 

242 help='Maximum traction vector length [Pa].') 

243 

244 def __init__(self, *args, **kwargs): 

245 super().__init__(*args, **kwargs) 

246 if self.rseed is None: 

247 self.rseed = num.random.randint(0, 2**31-1) 

248 self._data = None 

249 

250 def _get_data(self, nx, ny): 

251 if self._data is None: 

252 rstate = num.random.RandomState(self.rseed) 

253 self._data = rstate.rand(nx, ny) 

254 

255 return self._data 

256 

257 def get_tractions(self, nx, ny, patches=None): 

258 if patches is None: 

259 raise AttributeError( 

260 'Patches needs to be given for this traction field.') 

261 npatches = nx * ny 

262 dx = -patches[0].al1 + patches[0].al2 

263 dy = -patches[0].aw1 + patches[0].aw2 

264 

265 # Create random data and get spectrum and power spectrum 

266 data = self._get_data(nx, ny) 

267 spec = num.fft.fftshift(num.fft.fft2(data)) 

268 power_spec = (num.abs(spec)/spec.size)**2 

269 

270 # Get 0-centered wavenumbers (k_rad == 0.) is in the centre 

271 kx = num.fft.fftshift(num.fft.fftfreq(nx, d=dx)) 

272 ky = num.fft.fftshift(num.fft.fftfreq(ny, d=dy)) 

273 k_rad = num.sqrt(ky[:, num.newaxis]**2 + kx[num.newaxis, :]**2) 

274 

275 # Define wavenumber bins 

276 k_bins = num.arange(0, num.max(k_rad), num.max(k_rad)/10.) 

277 

278 # Set amplitudes within wavenumber bins to power_spec * 1 / k_max 

279 amps = num.zeros(k_rad.shape) 

280 amps[k_rad == 0.] = 1. 

281 

282 for i in range(k_bins.size-1): 

283 k_min = k_bins[i] 

284 k_max = k_bins[i+1] 

285 r = num.logical_and(k_rad > k_min, k_rad <= k_max) 

286 amps[r] = power_spec.T[r] 

287 amps = num.sqrt(amps * data.size * num.pi * 4) 

288 

289 amps[k_rad > k_bins.max()] = power_spec.ravel()[num.argmax(power_spec)] 

290 

291 # Multiply spectrum by amplitudes and inverse fft into demeaned noise 

292 spec *= amps.T 

293 

294 tractions = num.abs(num.fft.ifft2(spec)) 

295 tractions -= num.mean(tractions) 

296 tractions *= self.traction_max / num.abs(tractions).max() 

297 

298 t = num.zeros((npatches, 3)) 

299 t[:, 0] = num.cos(self.rake*d2r) * tractions.ravel(order='C') 

300 t[:, 1] = num.sin(self.rake*d2r) * tractions.ravel(order='C') 

301 

302 return t 

303 

304 

305class RectangularTaper(AbstractTractionField): 

306 width = Float.T( 

307 default=.2, 

308 help='Width of the taper as a fraction of the plane.') 

309 

310 type = StringChoice.T( 

311 choices=('tukey', ), 

312 default='tukey', 

313 help='Type of the taper, default: "tukey".') 

314 

315 def get_tractions(self, nx, ny, patches=None): 

316 if self.type == 'tukey': 

317 x = tukey_window(nx, self.width) 

318 y = tukey_window(ny, self.width) 

319 return (x[:, num.newaxis] * y).ravel()[:, num.newaxis] 

320 

321 raise AttributeError('Unknown type: %s' % self.type) 

322 

323 

324class DepthTaper(AbstractTractionField): 

325 depth_start = Float.T( 

326 help='Depth where the taper begins [m].') 

327 

328 depth_stop = Float.T( 

329 help='Depth where taper ends and drops to zero [m].') 

330 

331 type = StringChoice.T( 

332 choices=('linear', ), 

333 default='linear', 

334 help='Type of the taper, default: "linear".') 

335 

336 def get_tractions(self, nx, ny, patches): 

337 assert self.depth_stop > self.depth_start 

338 depths = num.array([p.depth for p in patches]) 

339 

340 if self.type == 'linear': 

341 slope = self.depth_stop - self.depth_start 

342 depths -= self.depth_stop 

343 depths /= -slope 

344 depths[depths > 1.] = 1. 

345 depths[depths < 0.] = 0. 

346 return depths[:, num.newaxis] 

347 

348 

349def plot_tractions(tractions, nx=15, ny=12, depth=10*km, component='strike'): 

350 ''' 

351 Plot traction model for quick inspection. 

352 

353 :param tractions: 

354 Traction field or traction composition to be displayed. 

355 :type tractions: 

356 :py:class:`pyrocko.gf.tractions.TractionField` 

357 

358 :param nx: 

359 Number of patches along strike. 

360 :type nx: 

361 int 

362 

363 :param ny: 

364 Number of patches down dip. 

365 :type ny: 

366 int 

367 

368 :param depth: 

369 Depth of the rupture plane center in [m]. 

370 :type depth: 

371 float 

372 

373 :param component: 

374 Choice of traction component to be shown. Available: ``'tx'`` (along 

375 strike), ``'ty'`` (up dip), ``'tz'`` (normal), ``'absolute'`` (vector 

376 length). 

377 :type component: 

378 str 

379 ''' 

380 import matplotlib.pyplot as plt 

381 from pyrocko.modelling.okada import OkadaSource 

382 

383 comp2idx = dict( 

384 tx=0, ty=1, tz=2) 

385 

386 source = OkadaSource( 

387 lat=0., 

388 lon=0., 

389 depth=depth, 

390 al1=-20*km, al2=20*km, 

391 aw1=-15*km, aw2=15*km, 

392 strike=120., dip=90., rake=90., 

393 slip=5.) 

394 

395 patches, _ = source.discretize(nx, ny) 

396 tractions = tractions.get_tractions(nx, ny, patches) 

397 

398 if component in comp2idx: 

399 tractions = tractions[:, comp2idx[component]].reshape(nx, ny) 

400 elif component == 'absolute': 

401 tractions = num.linalg.norm(tractions, axis=1).reshape(nx, ny) 

402 else: 

403 raise ValueError('Given component is not valid.') 

404 

405 fig = plt.figure() 

406 ax = fig.gca() 

407 

408 ax.imshow(tractions) 

409 

410 plt.show() 

411 

412 

413__all__ = [ 

414 'AbstractTractionField', 

415 'TractionField', 

416 'TractionComposition', 

417 'HomogeneousTractions', 

418 'DirectedTractions', 

419 'FractalTractions', 

420 'SelfSimilarTractions', 

421 'RectangularTaper', 

422 'DepthTaper', 

423 'plot_tractions']