Coverage for /usr/local/lib/python3.11/dist-packages/grond/problems/dynamic_rupture/problem.py: 30%

351 statements  

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

1import numpy as num 

2import logging 

3 

4from pyrocko import gf, util 

5from pyrocko.guts import String, Float, Dict, Int, Bool, StringChoice 

6from pyrocko.guts_array import Array 

7from pyrocko.gf.seismosizer import map_anchor 

8 

9from grond.meta import expand_template, Parameter, has_get_plot_classes, \ 

10 GrondError 

11 

12from ..base import Problem, ProblemConfig 

13from .. import CMTProblem 

14 

15guts_prefix = 'grond' 

16 

17logger = logging.getLogger('grond.problems.dynamic_rupture.problem') 

18km = 1e3 

19d2r = num.pi / 180. 

20as_km = dict(scale_factor=km, scale_unit='km') 

21as_gpa = dict(scale_factor=1e9, scale_unit='GPa') 

22 

23 

24class DynamicRuptureProblemConfig(ProblemConfig): 

25 

26 ranges = Dict.T(String.T(), gf.Range.T()) 

27 

28 anchor = StringChoice.T( 

29 choices=tuple(map_anchor.keys()), 

30 default='top') 

31 

32 decimation_factor = Int.T( 

33 default=1, 

34 help='Decimation factor of the discretized sub-faults') 

35 distance_min = Float.T(default=0.) 

36 adaptive_resolution = StringChoice.T( 

37 choices=('off', 'linear', 'uniqueness'), 

38 default='off') 

39 adaptive_start = Int.T( 

40 default=0) 

41 

42 pure_shear = Bool.T( 

43 default=True) 

44 smooth_rupture = Bool.T( 

45 default=False) 

46 

47 patch_mask = Array.T( 

48 optional=True, 

49 shape=(None,), 

50 dtype=bool, 

51 serialize_as='list') 

52 

53 point_source_target_balancing = Bool.T( 

54 default=False, 

55 help='If ``True``, target balancing (if used) is performed on a ' 

56 'moment tensor point source at the events location. It increases ' 

57 'the speed, but might lead to not fully optimized target weights.' 

58 ) 

59 

60 def get_problem(self, event, target_groups, targets): 

61 if self.decimation_factor != 1: 

62 logger.warn( 

63 'Decimation factor for dynamic rupture source set to %i.' 

64 ' Results may be inaccurate.' % self.decimation_factor) 

65 

66 base_source = gf.PseudoDynamicRupture.from_pyrocko_event( 

67 event, 

68 anchor=self.anchor, 

69 nx=int(self.ranges['nx'].start), 

70 ny=int(self.ranges['ny'].start), 

71 decimation_factor=self.decimation_factor, 

72 pure_shear=self.pure_shear, 

73 smooth_rupture=self.smooth_rupture, 

74 patch_mask=self.patch_mask, 

75 stf=None) 

76 

77 subs = dict( 

78 event_name=event.name, 

79 event_time=util.time_to_str(event.time)) 

80 

81 cmt_problem = None 

82 if self.point_source_target_balancing: 

83 base_source_cmt = gf.MTSource.from_pyrocko_event(event) 

84 

85 stf = gf.HalfSinusoidSTF() 

86 stf.duration = event.duration or 0.0 

87 

88 base_source_cmt.stf = stf 

89 

90 ranges = dict( 

91 time=self.ranges['time'], 

92 north_shift=self.ranges['north_shift'], 

93 east_shift=self.ranges['east_shift'], 

94 depth=self.ranges['depth'], 

95 magnitude=gf.Range( 

96 start=event.magnitude - 1., 

97 stop=event.magnitude + 1.), 

98 duration=gf.Range(start=0., stop=stf.duration * 2.), 

99 rmnn=gf.Range(start=-1.41421, stop=1.41421), 

100 rmee=gf.Range(start=-1.41421, stop=1.41421), 

101 rmdd=gf.Range(start=-1.41421, stop=1.41421), 

102 rmne=gf.Range(start=-1., stop=1.), 

103 rmnd=gf.Range(start=-1., stop=1.), 

104 rmed=gf.Range(start=-1., stop=1.)) 

105 

106 cmt_problem = CMTProblem( 

107 name=expand_template(self.name_template, subs), 

108 base_source=base_source_cmt, 

109 distance_min=self.distance_min, 

110 target_groups=target_groups, 

111 targets=targets, 

112 ranges=ranges, 

113 mt_type='dc', 

114 stf_type='HalfSinusoidSTF', 

115 norm_exponent=self.norm_exponent) 

116 

117 problem = DynamicRuptureProblem( 

118 name=expand_template(self.name_template, subs), 

119 base_source=base_source, 

120 distance_min=self.distance_min, 

121 target_groups=target_groups, 

122 targets=targets, 

123 ranges=self.ranges, 

124 norm_exponent=self.norm_exponent, 

125 adaptive_resolution=self.adaptive_resolution, 

126 adaptive_start=self.adaptive_start, 

127 cmt_problem=cmt_problem) 

128 

129 return problem 

130 

131 

132@has_get_plot_classes 

133class DynamicRuptureProblem(Problem): 

134 

135 problem_parameters = [ 

136 Parameter('east_shift', 'm', label='Easting', **as_km), 

137 Parameter('north_shift', 'm', label='Northing', **as_km), 

138 Parameter('depth', 'm', label='Depth', **as_km), 

139 Parameter('length', 'm', label='Length', **as_km), 

140 Parameter('width', 'm', label='Width', **as_km), 

141 Parameter('slip', 'm', label='Peak slip', optional=True), 

142 Parameter('strike', 'deg', label='Strike'), 

143 Parameter('dip', 'deg', label='Dip'), 

144 Parameter('rake', 'deg', label='Rake'), 

145 Parameter('gamma', 'vr/vs', label=r'$\gamma$'), 

146 Parameter('nx', 'patches', label='nx'), 

147 Parameter('ny', 'patches', label='ny') 

148 ] 

149 

150 problem_waveform_parameters = [ 

151 Parameter('nucleation_x', 'offset', label='Nucleation X'), 

152 Parameter('nucleation_y', 'offset', label='Nucleation Y'), 

153 Parameter('time', 's', label='Time'), 

154 ] 

155 

156 dependants = [ 

157 Parameter('traction', 'Pa', label='Maximum traction'), 

158 Parameter('magnitude', label='Magnitude'), 

159 ] 

160 

161 distance_min = Float.T(default=0.0) 

162 adaptive_resolution = String.T() 

163 adaptive_start = Int.T() 

164 

165 cmt_problem = Problem.T(optional=True) 

166 

167 def __init__(self, **kwargs): 

168 Problem.__init__(self, **kwargs) 

169 self.deps_cache = {} 

170 self.minmax_cache = {} 

171 

172 def set_engine(self, engine): 

173 self._engine = engine 

174 

175 if self.cmt_problem is not None: 

176 self.cmt_problem.set_engine(engine) 

177 

178 def pack(self, source): 

179 arr = self.get_parameter_array(source) 

180 for ip, p in enumerate(self.parameters): 

181 if p.name == 'time': 

182 arr[ip] -= self.base_source.time 

183 return arr 

184 

185 def get_source(self, x): 

186 src = self.base_source 

187 

188 d = self.get_parameter_dict(x) 

189 p = {k: float(self.ranges[k].make_relative(src[k], d[k])) 

190 for k in src.keys() if k in d} 

191 

192 p['nx'] = int(p['nx']) 

193 p['ny'] = int(p['ny']) 

194 

195 if p['nx'] != src.nx or p['ny'] != src.ny: 

196 logger.info('refining patches to %dx%d', p['nx'], p['ny']) 

197 

198 src.nx = p['nx'] 

199 src.ny = p['ny'] 

200 

201 return src.clone(**p) 

202 

203 def random_uniform(self, xbounds, rstate, fixed_magnitude=None): 

204 if fixed_magnitude is not None: 

205 raise GrondError( 

206 'Setting fixed magnitude in random model generation not ' 

207 'supported for this type of problem.') 

208 

209 x = num.zeros(self.nparameters) 

210 for i in range(self.nparameters): 

211 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1]) 

212 

213 return x 

214 

215 def preconstrain(self, x, optimiser=None): 

216 nx = self.ranges['nx'] 

217 ny = self.ranges['ny'] 

218 

219 if optimiser and self.adaptive_resolution == 'linear' and \ 

220 optimiser.iiter >= self.adaptive_start: 

221 

222 progress = (optimiser.iiter - self.adaptive_start - 1) / \ 

223 (optimiser.niterations - self.adaptive_start) 

224 

225 nx = num.floor( 

226 nx.start + progress * (nx.stop - nx.start + 1)) 

227 ny = num.floor( 

228 ny.start + progress * (ny.stop - ny.start + 1)) 

229 

230 elif optimiser and self.adaptive_resolution == 'uniqueness' and \ 

231 optimiser.iiter >= self.adaptive_start: 

232 raise NotImplementedError 

233 

234 else: 

235 nx = nx.start 

236 ny = ny.start 

237 

238 idx_nx = self.get_parameter_index('nx') 

239 idx_ny = self.get_parameter_index('ny') 

240 x[idx_nx] = int(round(nx)) 

241 x[idx_ny] = int(round(ny)) 

242 

243 return x 

244 

245 def make_dependant(self, xs, pname, store=None): 

246 cache = self.deps_cache 

247 if xs.ndim == 1: 

248 return self.make_dependant(xs[num.newaxis, :], pname)[0] 

249 

250 if pname not in self.dependant_names: 

251 raise KeyError(pname) 

252 

253 if store is None: 

254 store_ids = self.get_gf_store_ids() 

255 engine = self.get_engine() 

256 store = engine.get_store(store_ids[0]) 

257 

258 y = num.zeros(xs.shape[0]) 

259 for i, x in enumerate(xs): 

260 k = tuple(x.tolist()) 

261 if k not in cache: 

262 source = self.get_source(x) 

263 

264 if num.isnan(source.slip): 

265 source.slip = 0. 

266 

267 source.discretize_basesource(store) 

268 

269 tractions = source.get_scaled_tractions(store) 

270 traction = max(num.linalg.norm(tractions, axis=1)) 

271 

272 magnitude = source.get_magnitude(store=store) 

273 

274 cache[k] = dict( 

275 traction=traction, 

276 magnitude=magnitude) 

277 

278 res = cache[k] 

279 

280 y[i] = res[pname] 

281 

282 return y 

283 

284 def get_dependant_bounds(self): 

285 cache = self.minmax_cache 

286 

287 x_max_traction = self.get_min_model() 

288 

289 for i, p in enumerate(self.problem_parameters): 

290 if p.name != 'slip': 

291 continue 

292 

293 x_max_traction[i] = self.ranges[p.name].stop 

294 

295 xs = (self.get_min_model(), self.get_max_model(), x_max_traction) 

296 

297 if not any(tuple(x.tolist()) in cache for x in xs): 

298 store_ids = self.get_gf_store_ids() 

299 engine = self.get_engine() 

300 store = engine.get_store(store_ids[0]) 

301 

302 for x in xs: 

303 source = self.get_source(x) 

304 source.discretize_basesource(store) 

305 

306 tractions = source.get_scaled_tractions(store) 

307 traction = max(num.linalg.norm(tractions, axis=1)) 

308 

309 magnitude = source.get_magnitude(store=store) 

310 

311 cache[tuple(x.tolist())] = dict( 

312 traction=traction, 

313 magnitude=magnitude 

314 ) 

315 

316 traction = [cache[tuple(x.tolist())]['traction'] for x in xs] 

317 magnitude = [cache[tuple(x.tolist())]['magnitude'] for x in xs] 

318 

319 out = [ 

320 (min(traction), max(traction)), 

321 (min(magnitude), max(magnitude))] 

322 

323 return out 

324 

325 @classmethod 

326 def get_plot_classes(cls): 

327 from . import plot 

328 plots = super(DynamicRuptureProblem, cls).get_plot_classes() 

329 plots.extend([ 

330 plot.DynamicRuptureSlipMap, 

331 plot.DynamicRuptureSTF, 

332 plot.DynamicRuptureMap]) 

333 return plots 

334 

335 

336class DoublePDRProblemConfig(ProblemConfig): 

337 

338 ranges = Dict.T(String.T(), gf.Range.T()) 

339 

340 decimation_factor = Int.T( 

341 default=1, 

342 help='Decimation factor of the discretized sub-faults') 

343 distance_min = Float.T(default=0.) 

344 adaptive_resolution = StringChoice.T( 

345 choices=('off', 'linear', 'uniqueness'), 

346 default='off') 

347 adaptive_start = Int.T( 

348 default=0) 

349 

350 pure_shear = Bool.T( 

351 default=True) 

352 smooth_rupture = Bool.T( 

353 default=False) 

354 

355 point_source_target_balancing = Bool.T( 

356 default=False, 

357 help='If ``True``, target balancing (if used) is performed on a ' 

358 'moment tensor point source at the events location. It increases ' 

359 'the speed, but might lead to not fully optimized target weights.' 

360 ) 

361 

362 def get_problem(self, event, target_groups, targets): 

363 if self.decimation_factor != 1: 

364 logger.warn( 

365 'Decimation factor for dynamic rupture source set to %i.' 

366 ' Results may be inaccurate.' % self.decimation_factor) 

367 

368 base_source = gf.DoublePDR.from_pyrocko_event( 

369 event, 

370 nx1=int(self.ranges['nx1'].start), 

371 ny1=int(self.ranges['ny1'].start), 

372 nx2=int(self.ranges['nx2'].start), 

373 ny2=int(self.ranges['ny2'].start), 

374 decimation_factor=self.decimation_factor, 

375 pure_shear=self.pure_shear, 

376 smooth_rupture=self.smooth_rupture, 

377 stf=None) 

378 

379 subs = dict( 

380 event_name=event.name, 

381 event_time=util.time_to_str(event.time)) 

382 

383 cmt_problem = None 

384 if self.point_source_target_balancing: 

385 base_source_cmt = gf.MTSource.from_pyrocko_event(event) 

386 

387 stf = gf.HalfSinusoidSTF() 

388 stf.duration = event.duration or 0.0 

389 

390 base_source_cmt.stf = stf 

391 

392 ranges = dict( 

393 time=self.ranges['time'], 

394 north_shift=self.ranges['north_shift'], 

395 east_shift=self.ranges['east_shift'], 

396 depth=self.ranges['depth'], 

397 magnitude=gf.Range( 

398 start=event.magnitude - 1., 

399 stop=event.magnitude + 1.), 

400 duration=gf.Range(start=0., stop=stf.duration * 2.), 

401 rmnn=gf.Range(start=-1.41421, stop=1.41421), 

402 rmee=gf.Range(start=-1.41421, stop=1.41421), 

403 rmdd=gf.Range(start=-1.41421, stop=1.41421), 

404 rmne=gf.Range(start=-1., stop=1.), 

405 rmnd=gf.Range(start=-1., stop=1.), 

406 rmed=gf.Range(start=-1., stop=1.)) 

407 

408 cmt_problem = CMTProblem( 

409 name=expand_template(self.name_template, subs), 

410 base_source=base_source_cmt, 

411 distance_min=self.distance_min, 

412 target_groups=target_groups, 

413 targets=targets, 

414 ranges=ranges, 

415 mt_type='dc', 

416 stf_type='HalfSinusoidSTF', 

417 norm_exponent=self.norm_exponent) 

418 

419 problem = DoublePDRProblem( 

420 name=expand_template(self.name_template, subs), 

421 base_source=base_source, 

422 distance_min=self.distance_min, 

423 target_groups=target_groups, 

424 targets=targets, 

425 ranges=self.ranges, 

426 norm_exponent=self.norm_exponent, 

427 adaptive_resolution=self.adaptive_resolution, 

428 adaptive_start=self.adaptive_start, 

429 cmt_problem=cmt_problem) 

430 

431 return problem 

432 

433 

434@has_get_plot_classes 

435class DoublePDRProblem(Problem): 

436 

437 problem_parameters = [ 

438 Parameter('east_shift', 'm', label='Easting', **as_km), 

439 Parameter('north_shift', 'm', label='Northing', **as_km), 

440 Parameter('depth', 'm', label='Depth', **as_km), 

441 Parameter('length1', 'm', label='Length 1', **as_km), 

442 Parameter('width1', 'm', label='Width 1', **as_km), 

443 Parameter('slip1', 'm', label='Peak slip 1', optional=True), 

444 Parameter('strike1', 'deg', label='Strike 1'), 

445 Parameter('dip1', 'deg', label='Dip 1'), 

446 Parameter('rake1', 'deg', label='Rake 1'), 

447 Parameter('gamma1', 'vr/vs', label=r'$\gamma$ 1'), 

448 Parameter('nx1', 'patches', label='nx 1'), 

449 Parameter('ny1', 'patches', label='ny 1'), 

450 Parameter('length2', 'm', label='Length 2', **as_km), 

451 Parameter('width2', 'm', label='Width 2', **as_km), 

452 Parameter('slip2', 'm', label='Peak slip 2', optional=True), 

453 Parameter('strike2', 'deg', label='Strike 2'), 

454 Parameter('dip2', 'deg', label='Dip 2'), 

455 Parameter('rake2', 'deg', label='Rake 2'), 

456 Parameter('gamma2', 'vr/vs', label=r'$\gamma$ 2'), 

457 Parameter('nx2', 'patches', label='nx 2'), 

458 Parameter('ny2', 'patches', label='ny 2'), 

459 Parameter('delta_time', 's', label='$\\Delta$ Time'), 

460 Parameter('delta_depth', 'm', label='$\\Delta$ Depth'), 

461 Parameter('azimuth', 'deg', label='Azimuth'), 

462 Parameter('distance', 'm', label='Distance'), 

463 Parameter('mix', label='Mix') 

464 ] 

465 

466 problem_waveform_parameters = [ 

467 Parameter('nucleation_x1', 'offset', label='Nucleation X 1'), 

468 Parameter('nucleation_y1', 'offset', label='Nucleation Y 1'), 

469 Parameter('nucleation_x2', 'offset', label='Nucleation X 2'), 

470 Parameter('nucleation_y2', 'offset', label='Nucleation Y 2'), 

471 Parameter('time', 's', label='Time'), 

472 ] 

473 

474 dependants = [] 

475 

476 distance_min = Float.T(default=0.0) 

477 adaptive_resolution = String.T() 

478 adaptive_start = Int.T() 

479 

480 cmt_problem = Problem.T(optional=True) 

481 

482 def set_engine(self, engine): 

483 self._engine = engine 

484 

485 if self.cmt_problem is not None: 

486 self.cmt_problem.set_engine(engine) 

487 

488 def pack(self, source): 

489 arr = self.get_parameter_array(source) 

490 for ip, p in enumerate(self.parameters): 

491 if p.name == 'time': 

492 arr[ip] -= self.base_source.time 

493 return arr 

494 

495 def get_source(self, x): 

496 src = self.base_source 

497 

498 d = self.get_parameter_dict(x) 

499 p = {k: float(self.ranges[k].make_relative(src[k], d[k])) 

500 for k in src.keys() if k in d} 

501 

502 p['nx1'] = int(p['nx1']) 

503 p['ny1'] = int(p['ny1']) 

504 

505 if p['nx1'] != src.nx1 or p['ny1'] != src.ny1: 

506 logger.info('refining patches to %dx%d', p['nx1'], p['ny1']) 

507 

508 src.nx1 = p['nx1'] 

509 src.ny1 = p['ny1'] 

510 

511 p['nx2'] = int(p['nx2']) 

512 p['ny2'] = int(p['ny2']) 

513 

514 if p['nx2'] != src.nx2 or p['ny2'] != src.ny2: 

515 logger.info('refining patches to %dx%d', p['nx2'], p['ny2']) 

516 

517 src.nx2 = p['nx2'] 

518 src.ny2 = p['ny2'] 

519 

520 return src.clone(**p) 

521 

522 def random_uniform(self, xbounds, rstate, fixed_magnitude=None): 

523 if fixed_magnitude is not None: 

524 raise GrondError( 

525 'Setting fixed magnitude in random model generation not ' 

526 'supported for this type of problem.') 

527 

528 x = num.zeros(self.nparameters) 

529 for i in range(self.nparameters): 

530 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1]) 

531 

532 return x 

533 

534 def preconstrain(self, x, optimiser=None): 

535 nx1 = self.ranges['nx1'] 

536 ny1 = self.ranges['ny1'] 

537 nx2 = self.ranges['nx2'] 

538 ny2 = self.ranges['ny2'] 

539 

540 if optimiser and self.adaptive_resolution == 'linear' and \ 

541 optimiser.iiter >= self.adaptive_start: 

542 

543 progress = (optimiser.iiter - self.adaptive_start - 1) / \ 

544 (optimiser.niterations - self.adaptive_start) 

545 

546 nx1 = num.floor( 

547 nx1.start + progress * (nx1.stop - nx1.start + 1)) 

548 ny1 = num.floor( 

549 ny1.start + progress * (ny1.stop - ny1.start + 1)) 

550 nx2 = num.floor( 

551 nx2.start + progress * (nx2.stop - nx2.start + 1)) 

552 ny2 = num.floor( 

553 ny2.start + progress * (ny2.stop - ny2.start + 1)) 

554 

555 elif optimiser and self.adaptive_resolution == 'uniqueness' and \ 

556 optimiser.iiter >= self.adaptive_start: 

557 raise NotImplementedError 

558 

559 else: 

560 nx1 = nx1.start 

561 ny1 = ny1.start 

562 nx2 = nx2.start 

563 ny2 = ny2.start 

564 

565 idx_nx1 = self.get_parameter_index('nx1') 

566 idx_ny1 = self.get_parameter_index('ny1') 

567 idx_nx2 = self.get_parameter_index('nx2') 

568 idx_ny2 = self.get_parameter_index('ny2') 

569 x[idx_nx1] = int(round(nx1)) 

570 x[idx_ny1] = int(round(ny1)) 

571 x[idx_nx2] = int(round(nx2)) 

572 x[idx_ny2] = int(round(ny2)) 

573 

574 return x 

575 

576 @classmethod 

577 def get_plot_classes(cls): 

578 plots = super(DoublePDRProblem, cls).get_plot_classes() 

579 

580 return plots 

581 

582 

583class TriplePDRProblemConfig(ProblemConfig): 

584 

585 ranges = Dict.T(String.T(), gf.Range.T()) 

586 

587 decimation_factor = Int.T( 

588 default=1, 

589 help='Decimation factor of the discretized sub-faults') 

590 distance_min = Float.T(default=0.) 

591 adaptive_resolution = StringChoice.T( 

592 choices=('off', 'linear', 'uniqueness'), 

593 default='off') 

594 adaptive_start = Int.T( 

595 default=0) 

596 

597 pure_shear = Bool.T( 

598 default=True) 

599 smooth_rupture = Bool.T( 

600 default=False) 

601 

602 def get_problem(self, event, target_groups, targets): 

603 if self.decimation_factor != 1: 

604 logger.warn( 

605 'Decimation factor for dynamic rupture source set to %i.' 

606 ' Results may be inaccurate.' % self.decimation_factor) 

607 

608 base_source = gf.TriplePDR.from_pyrocko_event( 

609 event, 

610 nx1=int(self.ranges['nx1'].start), 

611 ny1=int(self.ranges['ny1'].start), 

612 nx2=int(self.ranges['nx2'].start), 

613 ny2=int(self.ranges['ny2'].start), 

614 nx3=int(self.ranges['nx3'].start), 

615 ny3=int(self.ranges['ny3'].start), 

616 decimation_factor=self.decimation_factor, 

617 pure_shear=self.pure_shear, 

618 smooth_rupture=self.smooth_rupture, 

619 stf=None) 

620 

621 subs = dict( 

622 event_name=event.name, 

623 event_time=util.time_to_str(event.time)) 

624 

625 cmt_problem = None 

626 

627 problem = TriplePDRProblem( 

628 name=expand_template(self.name_template, subs), 

629 base_source=base_source, 

630 distance_min=self.distance_min, 

631 target_groups=target_groups, 

632 targets=targets, 

633 ranges=self.ranges, 

634 norm_exponent=self.norm_exponent, 

635 adaptive_resolution=self.adaptive_resolution, 

636 adaptive_start=self.adaptive_start, 

637 cmt_problem=cmt_problem) 

638 

639 return problem 

640 

641 

642@has_get_plot_classes 

643class TriplePDRProblem(Problem): 

644 

645 problem_parameters = [ 

646 Parameter('delta_time1', 's', label=r'$\Delta$Time 1'), 

647 Parameter('north_shift1', 'm', label='Northing 1', **as_km), 

648 Parameter('east_shift1', 'm', label='Easting 1', **as_km), 

649 Parameter('depth1', 'm', label='Depth 1', **as_km), 

650 Parameter('length1', 'm', label='Length 1', **as_km), 

651 Parameter('width1', 'm', label='Width 1', **as_km), 

652 Parameter('slip1', 'm', label='Peak slip 1', optional=True), 

653 Parameter('strike1', 'deg', label='Strike 1'), 

654 Parameter('dip1', 'deg', label='Dip 1'), 

655 Parameter('rake1', 'deg', label='Rake 1'), 

656 Parameter('gamma1', 'vr/vs', label=r'$\gamma$ 1'), 

657 Parameter('nx1', 'patches', label='nx 1'), 

658 Parameter('ny1', 'patches', label='ny 1'), 

659 Parameter('delta_time2', 's', label=r'$\Delta$Time 2'), 

660 Parameter('north_shift2', 'm', label='Northing 2', **as_km), 

661 Parameter('east_shift2', 'm', label='Easting 2', **as_km), 

662 Parameter('depth2', 'm', label='Depth 2', **as_km), 

663 Parameter('length2', 'm', label='Length 2', **as_km), 

664 Parameter('width2', 'm', label='Width 2', **as_km), 

665 Parameter('slip2', 'm', label='Peak slip 2', optional=True), 

666 Parameter('strike2', 'deg', label='Strike 2'), 

667 Parameter('dip2', 'deg', label='Dip 2'), 

668 Parameter('rake2', 'deg', label='Rake 2'), 

669 Parameter('gamma2', 'vr/vs', label=r'$\gamma$ 2'), 

670 Parameter('nx2', 'patches', label='nx 2'), 

671 Parameter('ny2', 'patches', label='ny 2'), 

672 Parameter('delta_time3', 's', label=r'$\Delta$Time 3'), 

673 Parameter('north_shift3', 'm', label='Northing 3', **as_km), 

674 Parameter('east_shift3', 'm', label='Easting 3', **as_km), 

675 Parameter('depth3', 'm', label='Depth 3', **as_km), 

676 Parameter('length3', 'm', label='Length 3', **as_km), 

677 Parameter('width3', 'm', label='Width 3', **as_km), 

678 Parameter('slip3', 'm', label='Peak slip 3', optional=True), 

679 Parameter('strike3', 'deg', label='Strike 3'), 

680 Parameter('dip3', 'deg', label='Dip 3'), 

681 Parameter('rake3', 'deg', label='Rake 3'), 

682 Parameter('gamma3', 'vr/vs', label=r'$\gamma$ 3'), 

683 Parameter('nx3', 'patches', label='nx 3'), 

684 Parameter('ny3', 'patches', label='ny 3') 

685 ] 

686 

687 problem_waveform_parameters = [ 

688 Parameter('nucleation_x1', 'offset', label='Nucleation X 1'), 

689 Parameter('nucleation_y1', 'offset', label='Nucleation Y 1'), 

690 Parameter('nucleation_x2', 'offset', label='Nucleation X 2'), 

691 Parameter('nucleation_y2', 'offset', label='Nucleation Y 2'), 

692 Parameter('nucleation_x3', 'offset', label='Nucleation X 3'), 

693 Parameter('nucleation_y3', 'offset', label='Nucleation Y 3') 

694 ] 

695 

696 dependants = [] 

697 

698 distance_min = Float.T(default=0.0) 

699 adaptive_resolution = String.T() 

700 adaptive_start = Int.T() 

701 

702 cmt_problem = Problem.T(optional=True) 

703 

704 def set_engine(self, engine): 

705 self._engine = engine 

706 

707 if self.cmt_problem is not None: 

708 self.cmt_problem.set_engine(engine) 

709 

710 def pack(self, source): 

711 return self.get_parameter_array(source) 

712 # return arr 

713 

714 def get_source(self, x): 

715 src = self.base_source 

716 

717 d = self.get_parameter_dict(x) 

718 p = {k: float(self.ranges[k].make_relative(src[k], d[k])) 

719 for k in src.keys() if k in d} 

720 

721 p['nx1'] = int(p['nx1']) 

722 p['ny1'] = int(p['ny1']) 

723 

724 if p['nx1'] != src.nx1 or p['ny1'] != src.ny1: 

725 logger.info('refining patches to %dx%d', p['nx1'], p['ny1']) 

726 

727 src.nx1 = p['nx1'] 

728 src.ny1 = p['ny1'] 

729 

730 p['nx2'] = int(p['nx2']) 

731 p['ny2'] = int(p['ny2']) 

732 

733 if p['nx2'] != src.nx2 or p['ny2'] != src.ny2: 

734 logger.info('refining patches to %dx%d', p['nx2'], p['ny2']) 

735 

736 src.nx2 = p['nx2'] 

737 src.ny2 = p['ny2'] 

738 

739 p['nx3'] = int(p['nx3']) 

740 p['ny3'] = int(p['ny3']) 

741 

742 if p['nx3'] != src.nx3 or p['ny3'] != src.ny3: 

743 logger.info('refining patches to %dx%d', p['nx3'], p['ny3']) 

744 

745 src.nx3 = p['nx3'] 

746 src.ny3 = p['ny3'] 

747 

748 return src.clone(**p) 

749 

750 def random_uniform(self, xbounds, rstate, fixed_magnitude=None): 

751 if fixed_magnitude is not None: 

752 raise GrondError( 

753 'Setting fixed magnitude in random model generation not ' 

754 'supported for this type of problem.') 

755 

756 x = num.zeros(self.nparameters) 

757 for i in range(self.nparameters): 

758 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1]) 

759 

760 return x 

761 

762 def preconstrain(self, x, optimiser=None): 

763 nx1 = self.ranges['nx1'] 

764 ny1 = self.ranges['ny1'] 

765 nx2 = self.ranges['nx2'] 

766 ny2 = self.ranges['ny2'] 

767 nx3 = self.ranges['nx3'] 

768 ny3 = self.ranges['ny3'] 

769 

770 if optimiser and self.adaptive_resolution == 'linear' and \ 

771 optimiser.iiter >= self.adaptive_start: 

772 

773 progress = (optimiser.iiter - self.adaptive_start - 1) / \ 

774 (optimiser.niterations - self.adaptive_start) 

775 

776 nx1 = num.floor( 

777 nx1.start + progress * (nx1.stop - nx1.start + 1)) 

778 ny1 = num.floor( 

779 ny1.start + progress * (ny1.stop - ny1.start + 1)) 

780 nx2 = num.floor( 

781 nx2.start + progress * (nx2.stop - nx2.start + 1)) 

782 ny2 = num.floor( 

783 ny2.start + progress * (ny2.stop - ny2.start + 1)) 

784 nx3 = num.floor( 

785 nx3.start + progress * (nx3.stop - nx3.start + 1)) 

786 ny3 = num.floor( 

787 ny3.start + progress * (ny3.stop - ny3.start + 1)) 

788 

789 elif optimiser and self.adaptive_resolution == 'uniqueness' and \ 

790 optimiser.iiter >= self.adaptive_start: 

791 raise NotImplementedError 

792 

793 else: 

794 nx1 = nx1.start 

795 ny1 = ny1.start 

796 nx2 = nx2.start 

797 ny2 = ny2.start 

798 nx3 = nx3.start 

799 ny3 = ny3.start 

800 

801 idx_nx1 = self.get_parameter_index('nx1') 

802 idx_ny1 = self.get_parameter_index('ny1') 

803 idx_nx2 = self.get_parameter_index('nx2') 

804 idx_ny2 = self.get_parameter_index('ny2') 

805 idx_nx3 = self.get_parameter_index('nx3') 

806 idx_ny3 = self.get_parameter_index('ny3') 

807 x[idx_nx1] = int(round(nx1)) 

808 x[idx_ny1] = int(round(ny1)) 

809 x[idx_nx2] = int(round(nx2)) 

810 x[idx_ny2] = int(round(ny2)) 

811 x[idx_nx3] = int(round(nx3)) 

812 x[idx_ny3] = int(round(ny3)) 

813 

814 return x 

815 

816 @classmethod 

817 def get_plot_classes(cls): 

818 plots = super(TriplePDRProblem, cls).get_plot_classes() 

819 

820 return plots 

821 

822 

823__all__ = ''' 

824 DynamicRuptureProblem 

825 DynamicRuptureProblemConfig 

826 DoublePDRProblem 

827 DoublePDRProblemConfig 

828 TriplePDRProblem 

829 TriplePDRProblemConfig 

830'''.split()