Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/fomosto/psgrn_pscmp.py: 61%

594 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2024-03-07 11:54 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6import logging 

7import warnings 

8import os 

9from os.path import join as pjoin 

10import signal 

11import shutil 

12import copy 

13 

14import math 

15import numpy as num 

16 

17from tempfile import mkdtemp 

18from subprocess import Popen, PIPE 

19 

20from pyrocko.guts import Float, Int, Tuple, List, Object, String 

21from pyrocko.model import Location 

22from pyrocko import gf, util, trace, cake 

23 

24 

25# how to call the programs 

26program_bins = { 

27 'pscmp.2008a': 'fomosto_pscmp2008a', 

28 'psgrn.2008a': 'fomosto_psgrn2008a' 

29} 

30 

31psgrn_displ_names = ('uz', 'ur', 'ut') 

32psgrn_stress_names = ('szz', 'srr', 'stt', 'szr', 'srt', 'stz') 

33psgrn_tilt_names = ('tr', 'tt', 'rot') 

34psgrn_gravity_names = ('gd', 'gr') 

35psgrn_components = 'ep ss ds cl'.split() 

36 

37km = 1000. 

38day = 3600. * 24 

39 

40guts_prefix = 'pf' 

41logger = logging.getLogger('pyrocko.fomosto.psgrn_pscmp') 

42 

43 

44def have_backend(): 

45 for cmd in [[exe] for exe in program_bins.values()]: 

46 try: 

47 p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE) 

48 (stdout, stderr) = p.communicate() 

49 

50 except OSError: 

51 return False 

52 

53 return True 

54 

55 

56def nextpow2(i): 

57 return 2 ** int(math.ceil(math.log(i) / math.log(2.))) 

58 

59 

60def str_float_vals(vals): 

61 return ' '.join('%e' % val for val in vals) 

62 

63 

64def str_int_vals(vals): 

65 return ' '.join('%i' % val for val in vals) 

66 

67 

68def str_str_vals(vals): 

69 return ' '.join("'%s'" % val for val in vals) 

70 

71 

72def cake_model_to_config(mod): 

73 srows = [] 

74 for ir, row in enumerate(mod.to_scanlines(get_burgers=True)): 

75 depth, vp, vs, rho, qp, qs, eta1, eta2, alpha = row 

76 # replace qs with etas = 0. 

77 row = [depth / km, vp / km, vs / km, rho, eta1, eta2, alpha] 

78 srows.append('%i %15s' % (ir + 1, str_float_vals(row))) 

79 

80 return '\n'.join(srows), len(srows) 

81 

82 

83class PsGrnSpatialSampling(Object): 

84 

85 ''' 

86 Definition of spatial sampling for PSGRN. 

87 

88 Note: attributes in this class use non-standard units [km] to be consistent 

89 with PSGRN text file input. Please read the documentation carefully. 

90 ''' 

91 

92 n_steps = Int.T(default=10) 

93 start_distance = Float.T(default=0.) # start sampling [km] from source 

94 end_distance = Float.T(default=100.) # end 

95 

96 def string_for_config(self): 

97 return '%i %15e %15e' % (self.n_steps, self.start_distance, 

98 self.end_distance) 

99 

100 

101class PsGrnConfig(Object): 

102 

103 ''' 

104 Configuration for PSGRN. 

105 

106 Note: attributes in this class use non-standard units [km] and [days] to 

107 be consistent with PSGRN text file input. Please read the documentation 

108 carefully. 

109 ''' 

110 

111 version = String.T(default='2008a') 

112 

113 sampling_interval = Float.T( 

114 default=1.0, 

115 help='Exponential sampling 1.- equidistant, > 1. decreasing sampling' 

116 ' with distance') 

117 n_time_samples = Int.T( 

118 optional=True, 

119 help='Number of temporal GF samples up to max_time. Has to be equal' 

120 ' to a power of 2! If not, next power of 2 is taken.') 

121 max_time = Float.T( 

122 optional=True, 

123 help='Maximum time [days] after seismic event.') 

124 

125 gf_depth_spacing = Float.T( 

126 optional=True, 

127 help='Depth spacing [km] for the observation points. ' 

128 'If not defined depth spacing from the `StoreConfig` is used') 

129 gf_distance_spacing = Float.T( 

130 optional=True, 

131 help='Spatial spacing [km] for the observation points. ' 

132 'If not defined distance spacing from the `StoreConfig` is used') 

133 observation_depth = Float.T( 

134 default=0., 

135 help='Depth of observation points [km]') 

136 

137 def items(self): 

138 return dict(self.T.inamevals(self)) 

139 

140 

141class PsGrnConfigFull(PsGrnConfig): 

142 

143 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True) 

144 psgrn_outdir = String.T(default='psgrn_green/') 

145 

146 distance_grid = PsGrnSpatialSampling.T(PsGrnSpatialSampling.D()) 

147 depth_grid = PsGrnSpatialSampling.T(PsGrnSpatialSampling.D()) 

148 

149 sw_source_regime = Int.T(default=1) # 1-continental, 0-ocean 

150 sw_gravity = Int.T(default=0) 

151 

152 accuracy_wavenumber_integration = Float.T(default=0.025) 

153 

154 displ_filenames = Tuple.T(3, String.T(), default=psgrn_displ_names) 

155 stress_filenames = Tuple.T(6, String.T(), default=psgrn_stress_names) 

156 tilt_filenames = Tuple.T(3, String.T(), psgrn_tilt_names) 

157 gravity_filenames = Tuple.T(2, String.T(), psgrn_gravity_names) 

158 

159 @staticmethod 

160 def example(): 

161 conf = PsGrnConfigFull() 

162 conf.earthmodel_1d = cake.load_model().extract(depth_max=100*km) 

163 conf.psgrn_outdir = 'TEST_psgrn_functions/' 

164 return conf 

165 

166 def string_for_config(self): 

167 

168 assert self.earthmodel_1d is not None 

169 

170 d = self.__dict__.copy() 

171 

172 model_str, nlines = cake_model_to_config(self.earthmodel_1d) 

173 d['n_t2'] = nextpow2(self.n_time_samples) 

174 d['n_model_lines'] = nlines 

175 d['model_lines'] = model_str 

176 

177 d['str_psgrn_outdir'] = "'%s'" % './' 

178 

179 d['str_displ_filenames'] = str_str_vals(self.displ_filenames) 

180 d['str_stress_filenames'] = str_str_vals(self.stress_filenames) 

181 d['str_tilt_filenames'] = str_str_vals(self.tilt_filenames) 

182 d['str_gravity_filenames'] = str_str_vals(self.gravity_filenames) 

183 

184 d['str_distance_grid'] = self.distance_grid.string_for_config() 

185 d['str_depth_grid'] = self.depth_grid.string_for_config() 

186 

187 template = '''# autogenerated PSGRN input by psgrn.py 

188#============================================================================= 

189# This is input file of FORTRAN77 program "psgrn08a" for computing responses 

190# (Green's functions) of a multi-layered viscoelastic halfspace to point 

191# dislocation sources buried at different depths. All results will be stored in 

192# the given directory and provide the necessary data base for the program 

193# "pscmp07a" for computing time-dependent deformation, geoid and gravity changes 

194# induced by an earthquake with extended fault planes via linear superposition. 

195# For more details, please read the accompanying READ.ME file. 

196# 

197# written by Rongjiang Wang 

198# GeoForschungsZentrum Potsdam 

199# e-mail: wang@gfz-potsdam.de 

200# phone +49 331 2881209 

201# fax +49 331 2881204 

202# 

203# Last modified: Potsdam, Jan, 2008 

204# 

205################################################################# 

206## ## 

207## Cylindrical coordinates (Z positive downwards!) are used. ## 

208## ## 

209## If not specified otherwise, SI Unit System is used overall! ## 

210## ## 

211################################################################# 

212# 

213#------------------------------------------------------------------------------ 

214# 

215# PARAMETERS FOR SOURCE-OBSERVATION CONFIGURATIONS 

216# ================================================ 

217# 1. the uniform depth of the observation points [km], switch for oceanic (0) 

218# or continental(1) earthquakes; 

219# 2. number of (horizontal) observation distances (> 1 and <= nrmax defined in 

220# psgglob.h), start and end distances [km], ratio (>= 1.0) between max. and 

221# min. sampling interval (1.0 for equidistant sampling); 

222# 3. number of equidistant source depths (>= 1 and <= nzsmax defined in 

223# psgglob.h), start and end source depths [km]; 

224# 

225# r1,r2 = minimum and maximum horizontal source-observation 

226# distances (r2 > r1). 

227# zs1,zs2 = minimum and maximum source depths (zs2 >= zs1 > 0). 

228# 

229# Note that the same sampling rates dr_min and dzs will be used later by the 

230# program "pscmp07a" for discretizing the finite source planes to a 2D grid 

231# of point sources. 

232#------------------------------------------------------------------------------ 

233 %(observation_depth)e %(sw_source_regime)i 

234 %(str_distance_grid)s %(sampling_interval)e 

235 %(str_depth_grid)s 

236#------------------------------------------------------------------------------ 

237# 

238# PARAMETERS FOR TIME SAMPLING 

239# ============================ 

240# 1. number of time samples (<= ntmax def. in psgglob.h) and time window [days]. 

241# 

242# Note that nt (> 0) should be power of 2 (the fft-rule). If nt = 1, the 

243# coseismic (t = 0) changes will be computed; If nt = 2, the coseismic 

244# (t = 0) and steady-state (t -> infinity) changes will be computed; 

245# Otherwise, time series for the given time samples will be computed. 

246# 

247#------------------------------------------------------------------------------ 

248 %(n_t2)i %(max_time)f 

249#------------------------------------------------------------------------------ 

250# 

251# PARAMETERS FOR WAVENUMBER INTEGRATION 

252# ===================================== 

253# 1. relative accuracy of the wave-number integration (suggested: 0.1 - 0.01) 

254# 2. factor (> 0 and < 1) for including influence of earth's gravity on the 

255# deformation field (e.g. 0/1 = without / with 100percent gravity effect). 

256#------------------------------------------------------------------------------ 

257 %(accuracy_wavenumber_integration)e 

258 %(sw_gravity)i 

259#------------------------------------------------------------------------------ 

260# 

261# PARAMETERS FOR OUTPUT FILES 

262# =========================== 

263# 

264# 1. output directory 

265# 2. file names for 3 displacement components (uz, ur, ut) 

266# 3. file names for 6 stress components (szz, srr, stt, szr, srt, stz) 

267# 4. file names for radial and tangential tilt components (as measured by a 

268# borehole tiltmeter), rigid rotation of horizontal plane, geoid and gravity 

269# changes (tr, tt, rot, gd, gr) 

270# 

271# Note that all file or directory names should not be longer than 80 

272# characters. Directory and subdirectoy names must be separated and ended 

273# by / (unix) or \ (dos)! All file names should be given without extensions 

274# that will be appended automatically by ".ep" for the explosion (inflation) 

275# source, ".ss" for the strike-slip source, ".ds" for the dip-slip source, 

276# and ".cl" for the compensated linear vector dipole source) 

277# 

278#------------------------------------------------------------------------------ 

279 %(str_psgrn_outdir)s 

280 %(str_displ_filenames)s 

281 %(str_stress_filenames)s 

282 %(str_tilt_filenames)s %(str_gravity_filenames)s 

283#------------------------------------------------------------------------------ 

284# 

285# GLOBAL MODEL PARAMETERS 

286# ======================= 

287# 1. number of data lines of the layered model (<= lmax as defined in psgglob.h) 

288# 

289# The surface and the upper boundary of the half-space as well as the 

290# interfaces at which the viscoelastic parameters are continuous, are all 

291# defined by a single data line; All other interfaces, at which the 

292# viscoelastic parameters are discontinuous, are all defined by two data 

293# lines (upper-side and lower-side values). This input format could also be 

294# used for a graphic plot of the layered model. Layers which have different 

295# parameter values at top and bottom, will be treated as layers with a 

296# constant gradient, and will be discretised to a number of homogeneous 

297# sublayers. Errors due to the discretisation are limited within about 

298# 5percent (changeable, see psgglob.h). 

299# 

300# 2.... parameters of the multilayered model 

301# 

302# Burgers rheology (a Kelvin-Voigt body and a Maxwell body in series 

303# connection) for relaxation of shear modulus is implemented. No relaxation 

304# of compressional modulus is considered. 

305# 

306# eta1 = transient viscosity (dashpot of the Kelvin-Voigt body; <= 0 means 

307# infinity value) 

308# eta2 = steady-state viscosity (dashpot of the Maxwell body; <= 0 means 

309# infinity value) 

310# alpha = ratio between the effective and the unrelaxed shear modulus 

311# = mu1/(mu1+mu2) (> 0 and <= 1) 

312# 

313# Special cases: 

314# (1) Elastic: eta1 and eta2 <= 0 (i.e. infinity); alpha meaningless 

315# (2) Maxwell body: eta1 <= 0 (i.e. eta1 = infinity) 

316# or alpha = 1 (i.e. mu1 = infinity) 

317# (3) Standard-Linear-Solid: eta2 <= 0 (i.e. infinity) 

318#------------------------------------------------------------------------------ 

319 %(n_model_lines)i |int: no_model_lines; 

320#------------------------------------------------------------------------------ 

321# no depth[km] vp[km/s] vs[km/s] rho[kg/m^3] eta1[Pa*s] eta2[Pa*s] alpha 

322#------------------------------------------------------------------------------ 

323%(model_lines)s 

324#=======================end of input=========================================== 

325''' # noqa 

326 return (template % d).encode('ascii') 

327 

328 

329class PsGrnError(gf.store.StoreError): 

330 pass 

331 

332 

333def remove_if_exists(fn, force=False): 

334 if os.path.exists(fn): 

335 if force: 

336 os.remove(fn) 

337 else: 

338 raise gf.CannotCreate('file %s already exists' % fn) 

339 

340 

341class PsGrnRunner(object): 

342 ''' 

343 Wrapper object to execute the program fomosto_psgrn. 

344 ''' 

345 

346 def __init__(self, outdir): 

347 outdir = os.path.abspath(outdir) 

348 if not os.path.exists(outdir): 

349 os.mkdir(outdir) 

350 self.outdir = outdir 

351 self.config = None 

352 

353 def run(self, config, force=False): 

354 ''' 

355 Run the program with the specified configuration. 

356 

357 :param config: :py:class:`PsGrnConfigFull` 

358 :param force: boolean, set true to overwrite existing output 

359 ''' 

360 self.config = config 

361 

362 input_fn = pjoin(self.outdir, 'input') 

363 

364 remove_if_exists(input_fn, force=force) 

365 

366 with open(input_fn, 'wb') as f: 

367 input_str = config.string_for_config() 

368 logger.debug('===== begin psgrn input =====\n' 

369 '%s===== end psgrn input =====' % input_str.decode()) 

370 f.write(input_str) 

371 program = program_bins['psgrn.%s' % config.version] 

372 

373 old_wd = os.getcwd() 

374 

375 os.chdir(self.outdir) 

376 

377 interrupted = [] 

378 

379 def signal_handler(signum, frame): 

380 os.kill(proc.pid, signal.SIGTERM) 

381 interrupted.append(True) 

382 

383 original = signal.signal(signal.SIGINT, signal_handler) 

384 try: 

385 try: 

386 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE) 

387 

388 except OSError: 

389 os.chdir(old_wd) 

390 raise PsGrnError( 

391 '''could not start psgrn executable: "%s" 

392Available fomosto backends and download links to the modelling codes are listed 

393on 

394 

395 https://pyrocko.org/docs/current/apps/fomosto/backends.html 

396 

397''' % program) 

398 

399 (output_str, error_str) = proc.communicate(b'input\n') 

400 

401 finally: 

402 signal.signal(signal.SIGINT, original) 

403 

404 if interrupted: 

405 raise KeyboardInterrupt() 

406 

407 logger.debug('===== begin psgrn output =====\n' 

408 '%s===== end psgrn output =====' % output_str.decode()) 

409 

410 errmess = [] 

411 if proc.returncode != 0: 

412 errmess.append( 

413 'psgrn had a non-zero exit state: %i' % proc.returncode) 

414 

415 if error_str: 

416 logger.warning( 

417 'psgrn emitted something via stderr: \n\n%s' 

418 % error_str.decode()) 

419 # errmess.append('psgrn emitted something via stderr') 

420 

421 if output_str.lower().find(b'error') != -1: 

422 errmess.append("the string 'error' appeared in psgrn output") 

423 

424 if errmess: 

425 os.chdir(old_wd) 

426 raise PsGrnError(''' 

427===== begin psgrn input ===== 

428%s===== end psgrn input ===== 

429===== begin psgrn output ===== 

430%s===== end psgrn output ===== 

431===== begin psgrn error ===== 

432%s===== end psgrn error ===== 

433%s 

434psgrn has been invoked as "%s" 

435in the directory %s'''.lstrip() % ( 

436 input_str.decode(), output_str.decode(), error_str.decode(), 

437 '\n'.join(errmess), program, self.outdir)) 

438 

439 self.psgrn_output = output_str 

440 self.psgrn_error = error_str 

441 

442 os.chdir(old_wd) 

443 

444 

445pscmp_displ_names = ('un', 'ue', 'ud') 

446pscmp_stress_names = ('snn', 'see', 'sdd', 'sne', 'snd', 'sed') 

447pscmp_tilt_names = ('tn', 'te', 'rot') 

448pscmp_gravity_names = ('gd', 'gr') 

449pscmp_all_names = pscmp_displ_names + pscmp_stress_names + pscmp_tilt_names + \ 

450 pscmp_gravity_names 

451 

452pscmp_component_mapping = { 

453 'displ': (pscmp_displ_names, (2, 5)), 

454 'stress': (pscmp_stress_names, (5, 11)), 

455 'tilt': (pscmp_tilt_names, (11, 14)), 

456 'gravity': (pscmp_gravity_names, (14, 16)), 

457 'all': (pscmp_all_names, (2, 16)), 

458 } 

459 

460 

461def dsin(value): 

462 return num.sin(value * num.pi / 180.) 

463 

464 

465def dcos(value): 

466 return num.cos(value * num.pi / 180.) 

467 

468 

469def distributed_fault_patches_to_config(patches): 

470 ''' 

471 Input: List of PsCmpRectangularSource(s) 

472 ''' 

473 srows = [] 

474 for i, patch in enumerate(patches): 

475 srows.append('%i %s' % (i + 1, patch.string_for_config())) 

476 

477 return '\n'.join(srows), len(srows) 

478 

479 

480class PsCmpObservation(Object): 

481 pass 

482 

483 

484class PsCmpScatter(PsCmpObservation): 

485 ''' 

486 Scattered observation points. 

487 ''' 

488 lats = List.T(Float.T(), optional=True, default=[10.4, 10.5]) 

489 lons = List.T(Float.T(), optional=True, default=[12.3, 13.4]) 

490 

491 def string_for_config(self): 

492 srows = [] 

493 for lat, lon in zip(self.lats, self.lons): 

494 srows.append('(%15f, %15f)' % (lat, lon)) 

495 

496 self.sw = 0 

497 return ' %i' % (len(srows)), '\n'.join(srows) 

498 

499 

500class PsCmpProfile(PsCmpObservation): 

501 ''' 

502 Calculation along observation profile. 

503 ''' 

504 n_steps = Int.T(default=10) 

505 slat = Float.T( 

506 default=0., 

507 help='Profile start latitude') 

508 slon = Float.T( 

509 default=0., 

510 help='Profile start longitude') 

511 elat = Float.T( 

512 default=0., 

513 help='Profile end latitude') 

514 elon = Float.T( 

515 default=0., 

516 help='Profile end longitude') 

517 distances = List.T( 

518 Float.T(), 

519 optional=True, 

520 help='Distances [m] for each point on profile from start to end.') 

521 

522 def string_for_config(self): 

523 self.sw = 1 

524 

525 return ' %i' % (self.n_steps), \ 

526 ' ( %15f, %15f ), ( %15f, %15f )' % ( 

527 self.slat, self.slon, self.elat, self.elon) 

528 

529 

530class PsCmpArray(PsCmpObservation): 

531 ''' 

532 Calculation on a grid. 

533 ''' 

534 n_steps_lat = Int.T(default=10) 

535 n_steps_lon = Int.T(default=10) 

536 slat = Float.T( 

537 default=0., 

538 help='Array start latitude') 

539 slon = Float.T( 

540 default=0., 

541 help='Array start longitude') 

542 elat = Float.T( 

543 default=0., 

544 help='Array end latitude') 

545 elon = Float.T( 

546 default=0., 

547 help='Array end longitude') 

548 

549 def string_for_config(self): 

550 self.sw = 2 

551 

552 return ' %i %15f %15f ' % ( 

553 self.n_steps_lat, self.slat, self.elat), \ 

554 ' %i %15f %15f ' % ( 

555 self.n_steps_lon, self.slon, self.elon) 

556 

557 

558class PsCmpRectangularSource(Location, gf.seismosizer.Cloneable): 

559 ''' 

560 Rectangular source for the input geometry of the active fault. 

561 

562 The default shift of the origin (:gattr:`pos_s`, :gattr:`pos_d`) with 

563 respect to the reference coordinates 

564 (:gattr:`~pyrocko.model.location.Location.lat`, 

565 :gattr:`~pyrocko.model.location.Location.lon`) is zero, which implies that 

566 the reference is the center of the fault plane! The calculation point is 

567 always the center of the fault-plane. Setting :gattr:`pos_s` or 

568 :gattr:`pos_d` moves the fault point with respect to the origin along 

569 strike and dip direction, respectively. 

570 ''' 

571 

572 length = Float.T( 

573 default=6.0 * km, 

574 help='Fault length [m].') 

575 width = Float.T( 

576 default=5.0 * km, 

577 help='Fault width [m].') 

578 strike = Float.T( 

579 default=0.0, 

580 help='Fault strike, clockwise from north [deg].') 

581 dip = Float.T( 

582 default=90.0, 

583 help='Fault dip, downward from horizontal [deg].') 

584 rake = Float.T( 

585 default=0.0, 

586 help='Rake angle, counter-clock wise from right-horizontal [deg].') 

587 torigin = Float.T( 

588 default=0.0, 

589 help='Origin time [s].') 

590 slip = Float.T( 

591 optional=True, 

592 default=1.0, 

593 help='Slip amount [m]') 

594 pos_s = Float.T( 

595 optional=True, 

596 default=None, 

597 help='Origin offset along strike [m].') 

598 pos_d = Float.T( 

599 optional=True, 

600 default=None, 

601 help='Origin offset along dip [m].') 

602 opening = Float.T( 

603 default=0.0, 

604 help='Fault opening [m].') 

605 

606 def update(self, **kwargs): 

607 ''' 

608 Change some of the source models parameters. 

609 

610 Example:: 

611 

612 >>> from pyrocko import gf 

613 >>> s = gf.DCSource() 

614 >>> s.update(strike=66., dip=33.) 

615 >>> print(s) 

616 --- !pf.DCSource 

617 depth: 0.0 

618 time: 1970-01-01 00:00:00 

619 magnitude: 6.0 

620 strike: 66.0 

621 dip: 33.0 

622 rake: 0.0 

623 

624 ''' 

625 for (k, v) in kwargs.items(): 

626 self[k] = v 

627 

628 @property 

629 def dip_slip(self): 

630 return float(self.slip * dsin(self.rake) * (-1)) 

631 

632 @property 

633 def strike_slip(self): 

634 return float(self.slip * dcos(self.rake)) 

635 

636 def string_for_config(self): 

637 

638 if self.pos_s or self.pos_d is None: 

639 self.pos_s = 0. 

640 self.pos_d = 0. 

641 

642 tempd = copy.deepcopy(self.__dict__) 

643 tempd['dip_slip'] = self.dip_slip 

644 tempd['strike_slip'] = self.strike_slip 

645 tempd['effective_lat'] = self.effective_lat 

646 tempd['effective_lon'] = self.effective_lon 

647 tempd['depth'] /= km 

648 tempd['length'] /= km 

649 tempd['width'] /= km 

650 

651 return '%(effective_lat)15f %(effective_lon)15f %(depth)15f' \ 

652 '%(length)15f %(width)15f %(strike)15f' \ 

653 '%(dip)15f 1 1 %(torigin)15f \n %(pos_s)15f %(pos_d)15f ' \ 

654 '%(strike_slip)15f %(dip_slip)15f %(opening)15f' % tempd 

655 

656 

657MTIso = { 

658 'nn': dict(strike=90., dip=90., rake=0., slip=0., opening=1.), 

659 'ee': dict(strike=0., dip=90., rake=0., slip=0., opening=1.), 

660 'dd': dict(strike=0., dip=0., rake=-90., slip=0., opening=1.), 

661 } 

662 

663MTDev = { 

664 'ne': dict(strike=90., dip=90., rake=180., slip=1., opening=0.), 

665 'nd': dict(strike=180., dip=0., rake=0., slip=1., opening=0.), 

666 'ed': dict(strike=270., dip=0., rake=0., slip=1., opening=0.), 

667 } 

668 

669 

670def get_nullification_factor(mu, lame_lambda): 

671 ''' 

672 Factor that needs to be multiplied to 2 of the tensile sources to 

673 eliminate two of the isotropic components. 

674 ''' 

675 return -lame_lambda / (2. * mu + 2. * lame_lambda) 

676 

677 

678def get_trace_normalization_factor(mu, lame_lambda, nullification): 

679 ''' 

680 Factor to be multiplied to elementary GF trace to have unit displacement. 

681 ''' 

682 return 1. / ((2. * mu) + lame_lambda + (2. * lame_lambda * nullification)) 

683 

684 

685def get_iso_scaling_factors(mu, lame_lambda): 

686 nullf = get_nullification_factor(mu, lame_lambda) 

687 scale = get_trace_normalization_factor(mu, lame_lambda, nullf) 

688 return nullf, scale 

689 

690 

691class PsCmpTensileSF(Location, gf.seismosizer.Cloneable): 

692 ''' 

693 Compound dislocation of 3 perpendicular, rectangular sources to approximate 

694 an opening single force couple. NED coordinate system! 

695 

696 Warning: Mixed standard [m] / non-standard [days] units are used in this 

697 class. Please read the documentation carefully. 

698 ''' 

699 

700 length = Float.T( 

701 default=1.0 * km, 

702 help='Length of source rectangle [m].') 

703 width = Float.T( 

704 default=1.0 * km, 

705 help='Width of source rectangle [m].') 

706 torigin = Float.T( 

707 default=0.0, 

708 help='Origin time [days]') 

709 idx = String.T( 

710 default='nn', 

711 help='Axis index for opening direction; "nn", "ee", or "dd"') 

712 

713 def to_rfs(self, nullification): 

714 

715 cmpd = [] 

716 for comp, mt in MTIso.items(): 

717 params = copy.deepcopy(mt) 

718 

719 if comp != self.idx: 

720 params = copy.deepcopy(mt) 

721 params['opening'] *= nullification 

722 

723 kwargs = copy.deepcopy(self.__dict__) 

724 kwargs.update(params) 

725 kwargs.pop('idx') 

726 kwargs.pop('_latlon') 

727 cmpd.append(PsCmpRectangularSource(**kwargs)) 

728 

729 return cmpd 

730 

731 

732class PsCmpShearSF(Location, gf.seismosizer.Cloneable): 

733 ''' 

734 Shear fault source model component. 

735 

736 Warning: Mixed standard [m] / non-standard [days] units are used in this 

737 class. Please read the documentation carefully. 

738 ''' 

739 

740 length = Float.T( 

741 default=1.0 * km, 

742 help='Length of source rectangle [m].') 

743 width = Float.T( 

744 default=1.0 * km, 

745 help='Width of source rectangle [m].') 

746 torigin = Float.T( 

747 default=0.0, 

748 help='Origin time [days]') 

749 idx = String.T( 

750 default='ne', 

751 help='Axis index for opening direction; "ne", "nd", or "ed"') 

752 

753 def to_rfs(self): 

754 kwargs = copy.deepcopy(self.__dict__) 

755 kwargs.pop('idx') 

756 kwargs.pop('_latlon') 

757 kwargs.update(MTDev[self.idx]) 

758 return [PsCmpRectangularSource(**kwargs)] 

759 

760 

761class PsCmpMomentTensor(Location, gf.seismosizer.Cloneable): 

762 ''' 

763 Mapping of Moment Tensor components to rectangular faults. 

764 Only one component at a time valid! NED coordinate system! 

765 

766 Warning: Mixed standard [m] / non-standard [days] units are used in this 

767 class. Please read the documentation carefully. 

768 ''' 

769 length = Float.T( 

770 default=1.0 * km, 

771 help='Length of source rectangle [m].') 

772 width = Float.T( 

773 default=1.0 * km, 

774 help='Width of source rectangle [m].') 

775 torigin = Float.T( 

776 default=0.0, 

777 help='Origin time [days]') 

778 idx = String.T( 

779 default='nn', 

780 help='Axis index for MT component; ' 

781 '"nn", "ee", "dd", "ne", "nd", or "ed".') 

782 

783 def to_rfs(self, nullification=-0.25): 

784 kwargs = copy.deepcopy(self.__dict__) 

785 kwargs.pop('_latlon') 

786 

787 if self.idx in MTIso: 

788 tsf = PsCmpTensileSF(**kwargs) 

789 return tsf.to_rfs(nullification) 

790 

791 elif self.idx in MTDev: 

792 ssf = PsCmpShearSF(**kwargs) 

793 return ssf.to_rfs() 

794 

795 else: 

796 raise Exception('MT component not supported!') 

797 

798 

799class PsCmpCoulombStress(Object): 

800 pass 

801 

802 

803class PsCmpCoulombStressMasterFault(PsCmpCoulombStress): 

804 friction = Float.T(default=0.7) 

805 skempton_ratio = Float.T(default=0.0) 

806 master_fault_strike = Float.T(default=300.) 

807 master_fault_dip = Float.T(default=15.) 

808 master_fault_rake = Float.T(default=90.) 

809 sigma1 = Float.T(default=1.e6, help='[Pa]') 

810 sigma2 = Float.T(default=-1.e6, help='[Pa]') 

811 sigma3 = Float.T(default=0.0, help='[Pa]') 

812 

813 def string_for_config(self): 

814 return '%(friction)15e %(skempton_ratio)15e %(master_fault_strike)15f'\ 

815 '%(master_fault_dip)15f %(master_fault_rake)15f'\ 

816 '%(sigma1)15e %(sigma2)15e %(sigma3)15e' % self.__dict__ 

817 

818 

819class PsCmpSnapshots(Object): 

820 ''' 

821 Snapshot time series definition. 

822 

823 Note: attributes in this class use non-standard units [days] to be 

824 consistent with PSCMP text file input. Please read the documentation 

825 carefully. 

826 ''' 

827 

828 tmin = Float.T( 

829 default=0.0, 

830 help='Time [days] after source time to start temporal sample' 

831 ' snapshots.') 

832 tmax = Float.T( 

833 default=1.0, 

834 help='Time [days] after source time to end temporal sample f.') 

835 deltatdays = Float.T( 

836 default=1.0, 

837 help='Sample period [days].') 

838 

839 @property 

840 def times(self): 

841 nt = int(num.ceil((self.tmax - self.tmin) / self.deltatdays)) 

842 return num.linspace(self.tmin, self.tmax, nt).tolist() 

843 

844 @property 

845 def deltat(self): 

846 return self.deltatdays * 24 * 3600 

847 

848 

849class PsCmpConfig(Object): 

850 

851 version = String.T(default='2008a') 

852 # scatter, profile or array 

853 observation = PsCmpObservation.T(default=PsCmpScatter.D()) 

854 

855 rectangular_fault_size_factor = Float.T( 

856 default=1., 

857 help='The size of the rectangular faults in the compound MT' 

858 ' :py:class:`PsCmpMomentTensor` is dependend on the horizontal' 

859 ' spacing of the GF database. This factor is applied to the' 

860 ' relationship in i.e. fault length & width = f * dx.') 

861 

862 snapshots = PsCmpSnapshots.T( 

863 optional=True) 

864 

865 rectangular_source_patches = List.T(PsCmpRectangularSource.T()) 

866 

867 def items(self): 

868 return dict(self.T.inamevals(self)) 

869 

870 

871class PsCmpConfigFull(PsCmpConfig): 

872 

873 pscmp_outdir = String.T(default='./') 

874 psgrn_outdir = String.T(default='./psgrn_functions/') 

875 

876 los_vector = Tuple.T(3, Float.T(), optional=True) 

877 

878 sw_los_displacement = Int.T(default=0) 

879 sw_coulomb_stress = Int.T(default=0) 

880 coulomb_master_field = PsCmpCoulombStress.T( 

881 optional=True, 

882 default=PsCmpCoulombStressMasterFault.D()) 

883 

884 displ_sw_output_types = Tuple.T(3, Int.T(), default=(0, 0, 0)) 

885 stress_sw_output_types = Tuple.T(6, Int.T(), default=(0, 0, 0, 0, 0, 0)) 

886 tilt_sw_output_types = Tuple.T(3, Int.T(), default=(0, 0, 0)) 

887 gravity_sw_output_types = Tuple.T(2, Int.T(), default=(0, 0)) 

888 

889 indispl_filenames = Tuple.T(3, String.T(), default=psgrn_displ_names) 

890 instress_filenames = Tuple.T(6, String.T(), default=psgrn_stress_names) 

891 intilt_filenames = Tuple.T(3, String.T(), default=psgrn_tilt_names) 

892 ingravity_filenames = Tuple.T(2, String.T(), default=psgrn_gravity_names) 

893 

894 outdispl_filenames = Tuple.T(3, String.T(), default=pscmp_displ_names) 

895 outstress_filenames = Tuple.T(6, String.T(), default=pscmp_stress_names) 

896 outtilt_filenames = Tuple.T(3, String.T(), default=pscmp_tilt_names) 

897 outgravity_filenames = Tuple.T(2, String.T(), default=pscmp_gravity_names) 

898 

899 snapshot_basefilename = String.T(default='snapshot') 

900 

901 @classmethod 

902 def example(cls): 

903 conf = cls() 

904 conf.psgrn_outdir = 'TEST_psgrn_functions/' 

905 conf.pscmp_outdir = 'TEST_pscmp_output/' 

906 conf.rectangular_source_patches = [PsCmpRectangularSource( 

907 lat=10., lon=10., slip=2., 

908 width=5., length=10., 

909 strike=45, dip=30, rake=-90)] 

910 conf.observation = PsCmpArray( 

911 slat=9.5, elat=10.5, n_steps_lat=150, 

912 slon=9.5, elon=10.5, n_steps_lon=150) 

913 return conf 

914 

915 def get_output_filenames(self, rundir): 

916 return [pjoin(rundir, 

917 self.snapshot_basefilename + '_' + str(i + 1) + '.txt') 

918 for i in range(len(self.snapshots.times))] 

919 

920 def string_for_config(self): 

921 d = self.__dict__.copy() 

922 

923 patches_str, n_patches = distributed_fault_patches_to_config( 

924 self.rectangular_source_patches) 

925 

926 d['patches_str'] = patches_str 

927 d['n_patches'] = n_patches 

928 

929 str_npoints, str_observation = self.observation.string_for_config() 

930 d['str_npoints'] = str_npoints 

931 d['str_observation'] = str_observation 

932 d['sw_observation_type'] = self.observation.sw 

933 

934 if self.snapshots.times: 

935 str_times_dummy = [] 

936 for i, time in enumerate(self.snapshots.times): 

937 str_times_dummy.append("%f '%s_%i.txt'\n" % ( 

938 time, self.snapshot_basefilename, i + 1)) 

939 

940 str_times_dummy.append('#') 

941 d['str_times_snapshots'] = ''.join(str_times_dummy) 

942 d['n_snapshots'] = len(str_times_dummy) - 1 

943 else: 

944 d['str_times_snapshots'] = '# ' 

945 d['n_snapshots'] = 0 

946 

947 if self.sw_los_displacement == 1: 

948 d['str_los_vector'] = str_float_vals(self.los_vector) 

949 else: 

950 d['str_los_vector'] = '' 

951 

952 if self.sw_coulomb_stress == 1: 

953 d['str_coulomb_master_field'] = \ 

954 self.coulomb_master_field.string_for_config() 

955 else: 

956 d['str_coulomb_master_field'] = '' 

957 

958 d['str_psgrn_outdir'] = "'%s'" % self.psgrn_outdir 

959 d['str_pscmp_outdir'] = "'%s'" % './' 

960 

961 d['str_indispl_filenames'] = str_str_vals(self.indispl_filenames) 

962 d['str_instress_filenames'] = str_str_vals(self.instress_filenames) 

963 d['str_intilt_filenames'] = str_str_vals(self.intilt_filenames) 

964 d['str_ingravity_filenames'] = str_str_vals(self.ingravity_filenames) 

965 

966 d['str_outdispl_filenames'] = str_str_vals(self.outdispl_filenames) 

967 d['str_outstress_filenames'] = str_str_vals(self.outstress_filenames) 

968 d['str_outtilt_filenames'] = str_str_vals(self.outtilt_filenames) 

969 d['str_outgravity_filenames'] = str_str_vals(self.outgravity_filenames) 

970 

971 d['str_displ_sw_output_types'] = str_int_vals( 

972 self.displ_sw_output_types) 

973 d['str_stress_sw_output_types'] = str_int_vals( 

974 self.stress_sw_output_types) 

975 d['str_tilt_sw_output_types'] = str_int_vals( 

976 self.tilt_sw_output_types) 

977 d['str_gravity_sw_output_types'] = str_int_vals( 

978 self.gravity_sw_output_types) 

979 

980 template = '''# autogenerated PSCMP input by pscmp.py 

981#=============================================================================== 

982# This is input file of FORTRAN77 program "pscmp08" for modeling post-seismic 

983# deformation induced by earthquakes in multi-layered viscoelastic media using 

984# the Green's function approach. The earthquke source is represented by an 

985# arbitrary number of rectangular dislocation planes. For more details, please 

986# read the accompanying READ.ME file. 

987# 

988# written by Rongjiang Wang 

989# GeoForschungsZentrum Potsdam 

990# e-mail: wang@gfz-potsdam.de 

991# phone +49 331 2881209 

992# fax +49 331 2881204 

993# 

994# Last modified: Potsdam, July, 2008 

995# 

996# References: 

997# 

998# (1) Wang, R., F. Lorenzo-Martin and F. Roth (2003), Computation of deformation 

999# induced by earthquakes in a multi-layered elastic crust - FORTRAN programs 

1000# EDGRN/EDCMP, Computer and Geosciences, 29(2), 195-207. 

1001# (2) Wang, R., F. Lorenzo-Martin and F. Roth (2006), PSGRN/PSCMP - a new code for 

1002# calculating co- and post-seismic deformation, geoid and gravity changes 

1003# based on the viscoelastic-gravitational dislocation theory, Computers and 

1004# Geosciences, 32, 527-541. DOI:10.1016/j.cageo.2005.08.006. 

1005# (3) Wang, R. (2005), The dislocation theory: a consistent way for including the 

1006# gravity effect in (visco)elastic plane-earth models, Geophysical Journal 

1007# International, 161, 191-196. 

1008# 

1009################################################################# 

1010## ## 

1011## Green's functions should have been prepared with the ## 

1012## program "psgrn08" before the program "pscmp08" is started. ## 

1013## ## 

1014## For local Cartesian coordinate system, the Aki's convention ## 

1015## is used, that is, x is northward, y is eastward, and z is ## 

1016## downward. ## 

1017## ## 

1018## If not specified otherwise, SI Unit System is used overall! ## 

1019## ## 

1020################################################################# 

1021#=============================================================================== 

1022# OBSERVATION ARRAY 

1023# ================= 

1024# 1. selection for irregular observation positions (= 0) or a 1D observation 

1025# profile (= 1) or a rectangular 2D observation array (= 2): iposrec 

1026# 

1027# IF (iposrec = 0 for irregular observation positions) THEN 

1028# 

1029# 2. number of positions: nrec 

1030# 

1031# 3. coordinates of the observations: (lat(i),lon(i)), i=1,nrec 

1032# 

1033# ELSE IF (iposrec = 1 for regular 1D observation array) THEN 

1034# 

1035# 2. number of position samples of the profile: nrec 

1036# 

1037# 3. the start and end positions: (lat1,lon1), (lat2,lon2) 

1038# 

1039# ELSE IF (iposrec = 2 for rectanglular 2D observation array) THEN 

1040# 

1041# 2. number of x samples, start and end values: nxrec, xrec1, xrec2 

1042# 

1043# 3. number of y samples, start and end values: nyrec, yrec1, yrec2 

1044# 

1045# sequence of the positions in output data: lat(1),lon(1); ...; lat(nx),lon(1); 

1046# lat(1),lon(2); ...; lat(nx),lon(2); ...; lat(1),lon(ny); ...; lat(nx),lon(ny). 

1047# 

1048# Note that the total number of observation positions (nrec or nxrec*nyrec) 

1049# should be <= NRECMAX (see pecglob.h)! 

1050#=============================================================================== 

1051 %(sw_observation_type)i 

1052%(str_npoints)s 

1053%(str_observation)s 

1054#=============================================================================== 

1055# OUTPUTS 

1056# ======= 

1057# 

1058# 1. select output for los displacement (only for snapshots, see below), x, y, 

1059# and z-cosines to the INSAR orbit: insar (1/0 = yes/no), xlos, ylos, zlos 

1060# 

1061# if this option is selected (insar = 1), the snapshots will include additional 

1062# data: 

1063# LOS_Dsp = los displacement to the given satellite orbit. 

1064# 

1065# 2. select output for Coulomb stress changes (only for snapshots, see below): 

1066# icmb (1/0 = yes/no), friction, Skempton ratio, strike, dip, and rake angles 

1067# [deg] describing the uniform regional master fault mechanism, the uniform 

1068# regional principal stresses: sigma1, sigma2 and sigma3 [Pa] in arbitrary 

1069# order (the orietation of the pre-stress field will be derived by assuming 

1070# that the master fault is optimally oriented according to Coulomb failure 

1071# criterion) 

1072# 

1073# if this option is selected (icmb = 1), the snapshots will include additional 

1074# data: 

1075# CMB_Fix, Sig_Fix = Coulomb and normal stress changes on master fault; 

1076# CMB_Op1/2, Sig_Op1/2 = Coulomb and normal stress changes on the two optimally 

1077# oriented faults; 

1078# Str_Op1/2, Dip_Op1/2, Slp_Op1/2 = strike, dip and rake angles of the two 

1079# optimally oriented faults. 

1080# 

1081# Note: the 1. optimally orieted fault is the one closest to the master fault. 

1082# 

1083# 3. output directory in char format: outdir 

1084# 

1085# 4. select outputs for displacement components (1/0 = yes/no): itout(i), i=1-3 

1086# 

1087# 5. the file names in char format for the x, y, and z components: 

1088# toutfile(i), i=1-3 

1089# 

1090# 6. select outputs for stress components (1/0 = yes/no): itout(i), i=4-9 

1091# 

1092# 7. the file names in char format for the xx, yy, zz, xy, yz, and zx components: 

1093# toutfile(i), i=4-9 

1094# 

1095# 8. select outputs for vertical NS and EW tilt components, block rotation, geoid 

1096# and gravity changes (1/0 = yes/no): itout(i), i=10-14 

1097# 

1098# 9. the file names in char format for the NS tilt (positive if borehole top 

1099# tilts to north), EW tilt (positive if borehole top tilts to east), block 

1100# rotation (clockwise positive), geoid and gravity changes: toutfile(i), i=10-14 

1101# 

1102# Note that all above outputs are time series with the time window as same 

1103# as used for the Green's functions 

1104# 

1105#10. number of scenario outputs ("snapshots": spatial distribution of all above 

1106# observables at given time points; <= NSCENMAX (see pscglob.h): nsc 

1107# 

1108#11. the time [day], and file name (in char format) for the 1. snapshot; 

1109#12. the time [day], and file name (in char format) for the 2. snapshot; 

1110#13. ... 

1111# 

1112# Note that all file or directory names should not be longer than 80 

1113# characters. Directories must be ended by / (unix) or \ (dos)! 

1114#=============================================================================== 

1115 %(sw_los_displacement)i %(str_los_vector)s 

1116 %(sw_coulomb_stress)i %(str_coulomb_master_field)s 

1117 %(str_pscmp_outdir)s 

1118 %(str_displ_sw_output_types)s 

1119 %(str_outdispl_filenames)s 

1120 %(str_stress_sw_output_types)s 

1121 %(str_outstress_filenames)s 

1122 %(str_tilt_sw_output_types)s %(str_gravity_sw_output_types)s 

1123 %(str_outtilt_filenames)s %(str_outgravity_filenames)s 

1124 %(n_snapshots)i 

1125%(str_times_snapshots)s 

1126#=============================================================================== 

1127# 

1128# GREEN'S FUNCTION DATABASE 

1129# ========================= 

1130# 1. directory where the Green's functions are stored: grndir 

1131# 

1132# 2. file names (without extensions!) for the 13 Green's functions: 

1133# 3 displacement komponents (uz, ur, ut): green(i), i=1-3 

1134# 6 stress components (szz, srr, stt, szr, srt, stz): green(i), i=4-9 

1135# radial and tangential components measured by a borehole tiltmeter, 

1136# rigid rotation around z-axis, geoid and gravity changes (tr, tt, rot, gd, gr): 

1137# green(i), i=10-14 

1138# 

1139# Note that all file or directory names should not be longer than 80 

1140# characters. Directories must be ended by / (unix) or \ (dos)! The 

1141# extensions of the file names will be automatically considered. They 

1142# are ".ep", ".ss", ".ds" and ".cl" denoting the explosion (inflation) 

1143# strike-slip, the dip-slip and the compensated linear vector dipole 

1144# sources, respectively. 

1145# 

1146#=============================================================================== 

1147 %(str_psgrn_outdir)s 

1148 %(str_indispl_filenames)s 

1149 %(str_instress_filenames)s 

1150 %(str_intilt_filenames)s %(str_ingravity_filenames)s 

1151#=============================================================================== 

1152# RECTANGULAR SUBFAULTS 

1153# ===================== 

1154# 1. number of subfaults (<= NSMAX in pscglob.h): ns 

1155# 

1156# 2. parameters for the 1. rectangular subfault: geographic coordinates 

1157# (O_lat, O_lon) [deg] and O_depth [km] of the local reference point on 

1158# the present fault plane, length (along strike) [km] and width (along down 

1159# dip) [km], strike [deg], dip [deg], number of equi-size fault patches along 

1160# the strike (np_st) and along the dip (np_di) (total number of fault patches 

1161# = np_st x np_di), and the start time of the rupture; the following data 

1162# lines describe the slip distribution on the present sub-fault: 

1163# 

1164# pos_s[km] pos_d[km] slip_strike[m] slip_downdip[m] opening[m] 

1165# 

1166# where (pos_s,pos_d) defines the position of the center of each patch in 

1167# the local coordinate system with the origin at the reference point: 

1168# pos_s = distance along the length (positive in the strike direction) 

1169# pos_d = distance along the width (positive in the down-dip direction) 

1170# 

1171# 

1172# 3. ... for the 2. subfault ... 

1173# ... 

1174# N 

1175# / 

1176# /| strike 

1177# +------------------------ 

1178# |\ p . \ W 

1179# :-\ i . \ i 

1180# | \ l . \ d 

1181# :90 \ S . \ t 

1182# |-dip\ . \ h 

1183# : \. | rake \ 

1184# Z ------------------------- 

1185# L e n g t h 

1186# 

1187# Simulation of a Mogi source: 

1188# (1) Calculate deformation caused by three small openning plates (each 

1189# causes a third part of the volume of the point inflation) located 

1190# at the same depth as the Mogi source but oriented orthogonal to 

1191# each other. 

1192# (2) Multiply the results by 3(1-nu)/(1+nu), where nu is the Poisson 

1193# ratio at the source depth. 

1194# The multiplication factor is the ratio of the seismic moment (energy) of 

1195# the Mogi source to that of the plate openning with the same volume change. 

1196#=============================================================================== 

1197# n_faults 

1198#------------------------------------------------------------------------------- 

1199 %(n_patches)i 

1200#------------------------------------------------------------------------------- 

1201# n O_lat O_lon O_depth length width strike dip np_st np_di start_time 

1202# [-] [deg] [deg] [km] [km] [km] [deg] [deg] [-] [-] [day] 

1203# pos_s pos_d slp_stk slp_ddip open 

1204# [km] [km] [m] [m] [m] 

1205#------------------------------------------------------------------------------- 

1206%(patches_str)s 

1207#================================end of input=================================== 

1208''' # noqa 

1209 return (template % d).encode('ascii') 

1210 

1211 

1212class PsGrnPsCmpConfig(Object): 

1213 ''' 

1214 Combined config Object of PsGrn and PsCmp. 

1215 

1216 Note: attributes in this class use non-standard units [days] to be 

1217 consistent with PSCMP text file input. Please read the documentation 

1218 carefully. 

1219 ''' 

1220 tmin_days = Float.T( 

1221 default=0., 

1222 help='Min. time in days') 

1223 tmax_days = Float.T( 

1224 default=1., 

1225 help='Max. time in days') 

1226 gf_outdir = String.T(default='psgrn_functions') 

1227 

1228 psgrn_config = PsGrnConfig.T(default=PsGrnConfig.D()) 

1229 pscmp_config = PsCmpConfig.T(default=PsCmpConfig.D()) 

1230 

1231 

1232class PsCmpError(gf.store.StoreError): 

1233 pass 

1234 

1235 

1236class Interrupted(gf.store.StoreError): 

1237 def __str__(self): 

1238 return 'Interrupted.' 

1239 

1240 

1241class PsCmpRunner(object): 

1242 ''' 

1243 Wrapper object to execute the program fomosto_pscmp with the specified 

1244 configuration. 

1245 

1246 :param tmp: string, path to the temporary directy where calculation 

1247 results are stored 

1248 :param keep_tmp: boolean, if True the result directory is kept 

1249 ''' 

1250 def __init__(self, tmp=None, keep_tmp=False): 

1251 if tmp is not None: 

1252 tmp = os.path.abspath(tmp) 

1253 self.tempdir = mkdtemp(prefix='pscmprun-', dir=tmp) 

1254 self.keep_tmp = keep_tmp 

1255 self.config = None 

1256 

1257 def run(self, config): 

1258 ''' 

1259 Run the program! 

1260 

1261 :param config: :py:class:`PsCmpConfigFull` 

1262 ''' 

1263 self.config = config 

1264 

1265 input_fn = pjoin(self.tempdir, 'input') 

1266 

1267 with open(input_fn, 'wb') as f: 

1268 input_str = config.string_for_config() 

1269 

1270 logger.debug('===== begin pscmp input =====\n' 

1271 '%s===== end pscmp input =====' % input_str.decode()) 

1272 

1273 f.write(input_str) 

1274 

1275 program = program_bins['pscmp.%s' % config.version] 

1276 

1277 old_wd = os.getcwd() 

1278 os.chdir(self.tempdir) 

1279 

1280 interrupted = [] 

1281 

1282 def signal_handler(signum, frame): 

1283 os.kill(proc.pid, signal.SIGTERM) 

1284 interrupted.append(True) 

1285 

1286 original = signal.signal(signal.SIGINT, signal_handler) 

1287 try: 

1288 try: 

1289 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE, 

1290 close_fds=True) 

1291 

1292 except OSError as err: 

1293 os.chdir(old_wd) 

1294 logger.error('OS error: {0}'.format(err)) 

1295 raise PsCmpError( 

1296 '''could not start pscmp executable: "%s" 

1297Available fomosto backends and download links to the modelling codes are listed 

1298on 

1299 

1300 https://pyrocko.org/docs/current/apps/fomosto/backends.html 

1301 

1302''' % program) 

1303 

1304 (output_str, error_str) = proc.communicate(b'input\n') 

1305 

1306 finally: 

1307 signal.signal(signal.SIGINT, original) 

1308 

1309 if interrupted: 

1310 raise KeyboardInterrupt() 

1311 

1312 logger.debug('===== begin pscmp output =====\n' 

1313 '%s===== end pscmp output =====' % output_str.decode()) 

1314 

1315 errmsg = [] 

1316 if proc.returncode != 0: 

1317 errmsg.append( 

1318 'pscmp had a non-zero exit state: %i' % proc.returncode) 

1319 

1320 if error_str: 

1321 errmsg.append('pscmp emitted something via stderr') 

1322 

1323 if output_str.lower().find(b'error') != -1: 

1324 errmsg.append("the string 'error' appeared in pscmp output") 

1325 

1326 if errmsg: 

1327 self.keep_tmp = True 

1328 

1329 os.chdir(old_wd) 

1330 raise PsCmpError(''' 

1331===== begin pscmp input ===== 

1332{pscmp_input}===== end pscmp input ===== 

1333===== begin pscmp output ===== 

1334{pscmp_output}===== end pscmp output ===== 

1335===== begin pscmp error ===== 

1336{pscmp_error}===== end pscmp error ===== 

1337{error_messages} 

1338pscmp has been invoked as "{call}" 

1339in the directory {dir}'''.format( 

1340 pscmp_input=input_str, 

1341 pscmp_output=output_str, 

1342 pscmp_error=error_str, 

1343 error_messages='\n'.join(errmsg), 

1344 call=program, 

1345 dir=self.tempdir) 

1346 .lstrip()) 

1347 

1348 self.pscmp_output = output_str 

1349 self.pscmp_error = error_str 

1350 

1351 os.chdir(old_wd) 

1352 

1353 def get_results(self, component='displ'): 

1354 ''' 

1355 Get the resulting components from the stored directory. 

1356 

1357 .. note:: 

1358 

1359 The z-component is downward positive. 

1360 

1361 :param component: 

1362 The component to retrieve from the result directory. Choices: 

1363 ``"displ"`` -- displacement, n x 3 array, 

1364 ``"stress"`` -- stresses n x 6 array, 

1365 ``"tilt'`` -- tilts n x 3 array, 

1366 ``"gravity'`` -- gravity n x 2 array, 

1367 ``"all"`` -- all the above together. 

1368 :type param: 

1369 str 

1370 ''' 

1371 assert self.config.snapshots is not None 

1372 fns = self.config.get_output_filenames(self.tempdir) 

1373 

1374 output = [] 

1375 for fn in fns: 

1376 if not os.path.exists(fn): 

1377 continue 

1378 

1379 data = num.loadtxt(fn, skiprows=1, dtype=float) 

1380 

1381 try: 

1382 _, idxs = pscmp_component_mapping[component] 

1383 except KeyError: 

1384 raise Exception('component %s not supported! Either: %s' % ( 

1385 component, ', '.join( 

1386 '"%s"' % k for k in pscmp_component_mapping.keys()))) 

1387 

1388 istart, iend = idxs 

1389 

1390 output.append(data[:, istart:iend]) 

1391 

1392 return output 

1393 

1394 def get_traces(self, component='displ'): 

1395 ''' 

1396 Load snapshot data array and return specified components. 

1397 Transform array component and receiver wise to list of 

1398 :py:class:`pyrocko.trace.Trace` 

1399 ''' 

1400 

1401 distances = self.config.observation.distances 

1402 

1403 output_list = self.get_results(component=component) 

1404 deltat = self.config.snapshots.deltat 

1405 

1406 nrec, ncomp = output_list[0].shape 

1407 

1408 outarray = num.dstack([num.zeros([nrec, ncomp])] + output_list) 

1409 

1410 comps, _ = pscmp_component_mapping[component] 

1411 

1412 outtraces = [] 

1413 for row in range(nrec): 

1414 for col, comp in enumerate(comps): 

1415 tr = trace.Trace( 

1416 '', '%04i' % row, '', comp, 

1417 tmin=0.0, deltat=deltat, ydata=outarray[row, col, :], 

1418 meta=dict(distance=distances[row])) 

1419 outtraces.append(tr) 

1420 

1421 return outtraces 

1422 

1423 def __del__(self): 

1424 if self.tempdir: 

1425 if not self.keep_tmp: 

1426 shutil.rmtree(self.tempdir) 

1427 self.tempdir = None 

1428 else: 

1429 logger.warning( 

1430 'not removing temporary directory: %s' % self.tempdir) 

1431 

1432 

1433class PsGrnCmpGFBuilder(gf.builder.Builder): 

1434 nsteps = 2 

1435 

1436 def __init__(self, store_dir, step, shared, block_size=None, tmp=None, 

1437 force=False): 

1438 

1439 self.store = gf.store.Store(store_dir, 'w') 

1440 

1441 storeconf = self.store.config 

1442 

1443 dummy_lat = 10.0 

1444 dummy_lon = 10.0 

1445 

1446 depths = storeconf.coords[0] 

1447 lats = num.ones_like(depths) * dummy_lat 

1448 lons = num.ones_like(depths) * dummy_lon 

1449 points = num.vstack([lats, lons, depths]).T 

1450 

1451 self.shear_moduli = storeconf.get_shear_moduli( 

1452 lat=dummy_lat, 

1453 lon=dummy_lon, 

1454 points=points, 

1455 interpolation='multilinear') 

1456 

1457 self.lambda_moduli = storeconf.get_lambda_moduli( 

1458 lat=dummy_lat, 

1459 lon=dummy_lon, 

1460 points=points, 

1461 interpolation='multilinear') 

1462 

1463 if step == 0: 

1464 block_size = (1, storeconf.nsource_depths, storeconf.ndistances) 

1465 else: 

1466 if block_size is None: 

1467 block_size = (1, 1, storeconf.ndistances) 

1468 

1469 if len(storeconf.ns) == 2: 

1470 block_size = block_size[1:] 

1471 

1472 gf.builder.Builder.__init__( 

1473 self, storeconf, step, block_size=block_size, force=force) 

1474 

1475 baseconf = self.store.get_extra('psgrn_pscmp') 

1476 

1477 cg = PsGrnConfigFull(**baseconf.psgrn_config.items()) 

1478 if cg.n_time_samples is None or cg.max_time is None: 

1479 deltat_days = 1. / storeconf.sample_rate / day 

1480 

1481 cg.n_time_samples = int(baseconf.tmax_days // deltat_days) 

1482 cg.max_time = baseconf.tmax_days 

1483 else: 

1484 warnings.warn( 

1485 'PsGrnConfig is defining n_times_samples and max_time,' 

1486 ' this is replaced by PsGrnPsCmpConfig tmin and tmax.', 

1487 FutureWarning) 

1488 

1489 cc = PsCmpConfigFull(**baseconf.pscmp_config.items()) 

1490 if cc.snapshots is None: 

1491 deltat_days = 1. / storeconf.sample_rate / day 

1492 

1493 cc.snapshots = PsCmpSnapshots( 

1494 tmin=baseconf.tmin_days, 

1495 tmax=baseconf.tmax_days, 

1496 deltatdays=deltat_days) 

1497 else: 

1498 warnings.warn( 

1499 'PsCmpConfig is defining snapshots,' 

1500 ' this is replaced by PsGrnPsCmpConfig tmin and tmax.', 

1501 FutureWarning) 

1502 

1503 cg.earthmodel_1d = storeconf.earthmodel_1d 

1504 

1505 gf_outpath = os.path.join(store_dir, baseconf.gf_outdir) 

1506 

1507 cg.psgrn_outdir = gf_outpath + '/' 

1508 cc.psgrn_outdir = gf_outpath + '/' 

1509 

1510 util.ensuredir(gf_outpath) 

1511 

1512 self.cg = cg 

1513 self.cc = cc 

1514 

1515 def cleanup(self): 

1516 self.store.close() 

1517 

1518 def work_block(self, iblock): 

1519 

1520 if len(self.store.config.ns) == 2: 

1521 (sz, firstx), (sz, lastx), (ns, nx) = \ 

1522 self.get_block_extents(iblock) 

1523 mu = self.shear_moduli[iblock] 

1524 lame_lambda = self.lambda_moduli[iblock] 

1525 

1526 rz = self.store.config.receiver_depth 

1527 else: 

1528 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \ 

1529 self.get_block_extents(iblock) 

1530 

1531 fc = self.store.config 

1532 cc = copy.deepcopy(self.cc) 

1533 cg = copy.deepcopy(self.cg) 

1534 

1535 dx = fc.distance_delta 

1536 

1537 logger.info( 

1538 'Starting step %i / %i, block %i / %i' % 

1539 (self.step + 1, self.nsteps, iblock + 1, self.nblocks)) 

1540 

1541 if self.step == 0: 

1542 if cg.gf_depth_spacing is None: 

1543 gf_depth_spacing = fc.source_depth_delta 

1544 else: 

1545 gf_depth_spacing = cg.gf_depth_spacing * km 

1546 

1547 n_steps_depth = int((fc.source_depth_max - fc.source_depth_min) / 

1548 gf_depth_spacing) + 1 

1549 

1550 if cg.gf_distance_spacing is None: 

1551 gf_distance_spacing = fc.distance_delta 

1552 else: 

1553 gf_distance_spacing = cg.gf_distance_spacing * km 

1554 n_steps_distance = int((fc.distance_max - fc.distance_min) / 

1555 gf_distance_spacing) + 1 

1556 

1557 cg.depth_grid = PsGrnSpatialSampling( 

1558 n_steps=n_steps_depth, 

1559 start_distance=fc.source_depth_min / km, 

1560 end_distance=fc.source_depth_max / km) 

1561 

1562 cg.distance_grid = PsGrnSpatialSampling( 

1563 n_steps=n_steps_distance, 

1564 start_distance=fc.distance_min / km, 

1565 end_distance=fc.distance_max / km) 

1566 

1567 runner = PsGrnRunner(outdir=self.cg.psgrn_outdir) 

1568 runner.run(cg, force=self.force) 

1569 

1570 else: 

1571 distances = num.linspace( 

1572 firstx, firstx + (nx - 1) * dx, nx).tolist() 

1573 

1574 # fomosto sample rate in s, pscmp takes days 

1575 deltatdays = 1. / (fc.sample_rate * 24. * 3600.) 

1576 cc.snapshots = PsCmpSnapshots( 

1577 tmin=0., tmax=cg.max_time, deltatdays=deltatdays) 

1578 cc.observation = PsCmpProfile( 

1579 slat=0. - 0.001 * cake.m2d, 

1580 slon=0.0, 

1581 elat=0. + distances[-1] * cake.m2d, 

1582 elon=0.0, 

1583 n_steps=len(distances), 

1584 distances=distances) 

1585 

1586 runner = PsCmpRunner(keep_tmp=False) 

1587 

1588 mtsize = float(dx * cc.rectangular_fault_size_factor) 

1589 

1590 Ai = 1. / num.power(mtsize, 2) 

1591 nullf, sf = get_iso_scaling_factors(mu=mu, lame_lambda=lame_lambda) 

1592 

1593 mui = 1. / mu 

1594 

1595 gfmapping = [ 

1596 (('nn',), 

1597 {'un': (0, Ai * sf), 'ud': (5, Ai * sf)}), 

1598 (('ne',), 

1599 {'ue': (3, 1 * Ai * mui)}), 

1600 (('nd',), 

1601 {'un': (1, 1 * Ai * mui), 'ud': (6, 1 * Ai * mui)}), 

1602 (('ed',), 

1603 {'ue': (4, 1 * Ai * mui)}), 

1604 (('dd',), 

1605 {'un': (2, Ai * sf), 'ud': (7, Ai * sf)}), 

1606 (('ee',), 

1607 {'un': (8, Ai * sf), 'ud': (9, Ai * sf)}), 

1608 ] 

1609 

1610 for mt, gfmap in gfmapping: 

1611 cc.rectangular_source_patches = [] 

1612 for idx in mt: 

1613 pmt = PsCmpMomentTensor( 

1614 lat=0. + 0.001 * dx * cake.m2d, 

1615 lon=0.0, 

1616 depth=float(sz), 

1617 width=mtsize, 

1618 length=mtsize, 

1619 idx=idx) 

1620 

1621 cc.rectangular_source_patches.extend( 

1622 pmt.to_rfs(nullf)) 

1623 

1624 runner.run(cc) 

1625 

1626 rawtraces = runner.get_traces() 

1627 

1628 interrupted = [] 

1629 

1630 def signal_handler(signum, frame): 

1631 interrupted.append(True) 

1632 

1633 original = signal.signal(signal.SIGINT, signal_handler) 

1634 self.store.lock() 

1635 duplicate_inserts = 0 

1636 

1637 try: 

1638 for itr, tr in enumerate(rawtraces): 

1639 if tr.channel in gfmap: 

1640 x = tr.meta['distance'] 

1641 ig, factor = gfmap[tr.channel] 

1642 

1643 if len(self.store.config.ns) == 2: 

1644 args = (sz, x, ig) 

1645 else: 

1646 args = (rz, sz, x, ig) 

1647 

1648 gf_tr = gf.store.GFTrace.from_trace(tr) 

1649 

1650 gf_tr.data *= factor 

1651 

1652 try: 

1653 self.store.put(args, gf_tr) 

1654 except gf.store.DuplicateInsert: 

1655 duplicate_inserts += 1 

1656 

1657 finally: 

1658 if duplicate_inserts: 

1659 logger.warning( 

1660 '%i insertions skipped (duplicates)' 

1661 % duplicate_inserts) 

1662 

1663 self.store.unlock() 

1664 signal.signal(signal.SIGINT, original) 

1665 

1666 if interrupted: 

1667 raise KeyboardInterrupt() 

1668 

1669 logger.info( 

1670 'Done with step %i / %i, block %i / %i' % ( 

1671 self.step + 1, self.nsteps, iblock + 1, self.nblocks)) 

1672 

1673 

1674def init(store_dir, variant): 

1675 if variant is None: 

1676 variant = '2008a' 

1677 

1678 if ('pscmp.' + variant) not in program_bins: 

1679 raise gf.store.StoreError('unsupported pscmp variant: %s' % variant) 

1680 

1681 if ('psgrn.' + variant) not in program_bins: 

1682 raise gf.store.StoreError('unsupported psgrn variant: %s' % variant) 

1683 

1684 c = PsGrnPsCmpConfig() 

1685 

1686 store_id = os.path.basename(os.path.realpath(store_dir)) 

1687 

1688 # Initialising a viscous mantle 

1689 cake_mod = cake.load_model(fn=None, crust2_profile=(54., 23.)) 

1690 mantle = cake_mod.material(z=45*km) 

1691 mantle.burger_eta1 = 5e17 

1692 mantle.burger_eta2 = 1e19 

1693 mantle.burger_alpha = 1. 

1694 

1695 config = gf.meta.ConfigTypeA( 

1696 id=store_id, 

1697 ncomponents=10, 

1698 sample_rate=1. / (3600. * 24.), 

1699 receiver_depth=0. * km, 

1700 source_depth_min=0. * km, 

1701 source_depth_max=15. * km, 

1702 source_depth_delta=.5 * km, 

1703 distance_min=0. * km, 

1704 distance_max=50. * km, 

1705 distance_delta=1. * km, 

1706 earthmodel_1d=cake_mod, 

1707 modelling_code_id='psgrn_pscmp.%s' % variant, 

1708 tabulated_phases=[]) # dummy list 

1709 

1710 c.validate() 

1711 config.validate() 

1712 return gf.store.Store.create_editables( 

1713 store_dir, 

1714 config=config, 

1715 extra={'psgrn_pscmp': c}) 

1716 

1717 

1718def build(store_dir, 

1719 force=False, 

1720 nworkers=None, 

1721 continue_=False, 

1722 step=None, 

1723 iblock=None): 

1724 

1725 return PsGrnCmpGFBuilder.build( 

1726 store_dir, force=force, nworkers=nworkers, continue_=continue_, 

1727 step=step, iblock=iblock)