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 2023-10-06 15:01 +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 Input parameters have to be in: 

563 [deg] for reference point (lat, lon) and angles (rake, strike, dip) 

564 [m] shifting with respect to reference position 

565 [m] for fault dimensions and source depth. The default shift of the 

566 origin (:py:attr`pos_s`, :py:attr:`pos_d`) with respect to the reference 

567 coordinates 

568 (lat, lon) is zero, which implies that the reference is the center of 

569 the fault plane! 

570 The calculation point is always the center of the fault-plane! 

571 Setting :py:attr`pos_s` or :py:attr`pos_d` moves the fault point with 

572 respect to the origin along strike and dip direction, respectively! 

573 ''' 

574 length = Float.T(default=6.0 * km) 

575 width = Float.T(default=5.0 * km) 

576 strike = Float.T(default=0.0) 

577 dip = Float.T(default=90.0) 

578 rake = Float.T(default=0.0) 

579 torigin = Float.T(default=0.0) 

580 

581 slip = Float.T(optional=True, default=1.0) 

582 

583 pos_s = Float.T(optional=True, default=None) 

584 pos_d = Float.T(optional=True, default=None) 

585 opening = Float.T(default=0.0) 

586 

587 def update(self, **kwargs): 

588 ''' 

589 Change some of the source models parameters. 

590 

591 Example:: 

592 

593 >>> from pyrocko import gf 

594 >>> s = gf.DCSource() 

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

596 >>> print(s) 

597 --- !pf.DCSource 

598 depth: 0.0 

599 time: 1970-01-01 00:00:00 

600 magnitude: 6.0 

601 strike: 66.0 

602 dip: 33.0 

603 rake: 0.0 

604 

605 ''' 

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

607 self[k] = v 

608 

609 @property 

610 def dip_slip(self): 

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

612 

613 @property 

614 def strike_slip(self): 

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

616 

617 def string_for_config(self): 

618 

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

620 self.pos_s = 0. 

621 self.pos_d = 0. 

622 

623 tempd = copy.deepcopy(self.__dict__) 

624 tempd['dip_slip'] = self.dip_slip 

625 tempd['strike_slip'] = self.strike_slip 

626 tempd['effective_lat'] = self.effective_lat 

627 tempd['effective_lon'] = self.effective_lon 

628 tempd['depth'] /= km 

629 tempd['length'] /= km 

630 tempd['width'] /= km 

631 

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

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

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

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

636 

637 

638MTIso = { 

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

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

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

642 } 

643 

644MTDev = { 

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

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

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

648 } 

649 

650 

651def get_nullification_factor(mu, lame_lambda): 

652 ''' 

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

654 eliminate two of the isotropic components. 

655 ''' 

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

657 

658 

659def get_trace_normalization_factor(mu, lame_lambda, nullification): 

660 ''' 

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

662 ''' 

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

664 

665 

666def get_iso_scaling_factors(mu, lame_lambda): 

667 nullf = get_nullification_factor(mu, lame_lambda) 

668 scale = get_trace_normalization_factor(mu, lame_lambda, nullf) 

669 return nullf, scale 

670 

671 

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

673 ''' 

674 Compound dislocation of 3 perpendicular, rectangular sources to approximate 

675 an opening single force couple. NED coordinate system! 

676 

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

678 class. Please read the documentation carefully. 

679 ''' 

680 

681 length = Float.T( 

682 default=1.0 * km, 

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

684 width = Float.T( 

685 default=1.0 * km, 

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

687 torigin = Float.T( 

688 default=0.0, 

689 help='Origin time [days]') 

690 idx = String.T( 

691 default='nn', 

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

693 

694 def to_rfs(self, nullification): 

695 

696 cmpd = [] 

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

698 params = copy.deepcopy(mt) 

699 

700 if comp != self.idx: 

701 params = copy.deepcopy(mt) 

702 params['opening'] *= nullification 

703 

704 kwargs = copy.deepcopy(self.__dict__) 

705 kwargs.update(params) 

706 kwargs.pop('idx') 

707 kwargs.pop('_latlon') 

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

709 

710 return cmpd 

711 

712 

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

714 ''' 

715 Shear fault source model component. 

716 

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

718 class. Please read the documentation carefully. 

719 ''' 

720 

721 length = Float.T( 

722 default=1.0 * km, 

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

724 width = Float.T( 

725 default=1.0 * km, 

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

727 torigin = Float.T( 

728 default=0.0, 

729 help='Origin time [days]') 

730 idx = String.T( 

731 default='ne', 

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

733 

734 def to_rfs(self): 

735 kwargs = copy.deepcopy(self.__dict__) 

736 kwargs.pop('idx') 

737 kwargs.pop('_latlon') 

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

739 return [PsCmpRectangularSource(**kwargs)] 

740 

741 

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

743 ''' 

744 Mapping of Moment Tensor components to rectangular faults. 

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

746 

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

748 class. Please read the documentation carefully. 

749 ''' 

750 length = Float.T( 

751 default=1.0 * km, 

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

753 width = Float.T( 

754 default=1.0 * km, 

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

756 torigin = Float.T( 

757 default=0.0, 

758 help='Origin time [days]') 

759 idx = String.T( 

760 default='nn', 

761 help='Axis index for MT component; ' 

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

763 

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

765 kwargs = copy.deepcopy(self.__dict__) 

766 kwargs.pop('_latlon') 

767 

768 if self.idx in MTIso: 

769 tsf = PsCmpTensileSF(**kwargs) 

770 return tsf.to_rfs(nullification) 

771 

772 elif self.idx in MTDev: 

773 ssf = PsCmpShearSF(**kwargs) 

774 return ssf.to_rfs() 

775 

776 else: 

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

778 

779 

780class PsCmpCoulombStress(Object): 

781 pass 

782 

783 

784class PsCmpCoulombStressMasterFault(PsCmpCoulombStress): 

785 friction = Float.T(default=0.7) 

786 skempton_ratio = Float.T(default=0.0) 

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

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

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

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

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

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

793 

794 def string_for_config(self): 

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

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

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

798 

799 

800class PsCmpSnapshots(Object): 

801 ''' 

802 Snapshot time series definition. 

803 

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

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

806 carefully. 

807 ''' 

808 

809 tmin = Float.T( 

810 default=0.0, 

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

812 ' snapshots.') 

813 tmax = Float.T( 

814 default=1.0, 

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

816 deltatdays = Float.T( 

817 default=1.0, 

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

819 

820 @property 

821 def times(self): 

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

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

824 

825 @property 

826 def deltat(self): 

827 return self.deltatdays * 24 * 3600 

828 

829 

830class PsCmpConfig(Object): 

831 

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

833 # scatter, profile or array 

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

835 

836 rectangular_fault_size_factor = Float.T( 

837 default=1., 

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

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

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

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

842 

843 snapshots = PsCmpSnapshots.T( 

844 optional=True) 

845 

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

847 

848 def items(self): 

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

850 

851 

852class PsCmpConfigFull(PsCmpConfig): 

853 

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

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

856 

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

858 

859 sw_los_displacement = Int.T(default=0) 

860 sw_coulomb_stress = Int.T(default=0) 

861 coulomb_master_field = PsCmpCoulombStress.T( 

862 optional=True, 

863 default=PsCmpCoulombStressMasterFault.D()) 

864 

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

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

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

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

869 

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

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

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

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

874 

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

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

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

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

879 

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

881 

882 @classmethod 

883 def example(cls): 

884 conf = cls() 

885 conf.psgrn_outdir = 'TEST_psgrn_functions/' 

886 conf.pscmp_outdir = 'TEST_pscmp_output/' 

887 conf.rectangular_source_patches = [PsCmpRectangularSource( 

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

889 width=5., length=10., 

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

891 conf.observation = PsCmpArray( 

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

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

894 return conf 

895 

896 def get_output_filenames(self, rundir): 

897 return [pjoin(rundir, 

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

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

900 

901 def string_for_config(self): 

902 d = self.__dict__.copy() 

903 

904 patches_str, n_patches = distributed_fault_patches_to_config( 

905 self.rectangular_source_patches) 

906 

907 d['patches_str'] = patches_str 

908 d['n_patches'] = n_patches 

909 

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

911 d['str_npoints'] = str_npoints 

912 d['str_observation'] = str_observation 

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

914 

915 if self.snapshots.times: 

916 str_times_dummy = [] 

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

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

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

920 

921 str_times_dummy.append('#') 

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

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

924 else: 

925 d['str_times_snapshots'] = '# ' 

926 d['n_snapshots'] = 0 

927 

928 if self.sw_los_displacement == 1: 

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

930 else: 

931 d['str_los_vector'] = '' 

932 

933 if self.sw_coulomb_stress == 1: 

934 d['str_coulomb_master_field'] = \ 

935 self.coulomb_master_field.string_for_config() 

936 else: 

937 d['str_coulomb_master_field'] = '' 

938 

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

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

941 

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

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

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

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

946 

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

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

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

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

951 

952 d['str_displ_sw_output_types'] = str_int_vals( 

953 self.displ_sw_output_types) 

954 d['str_stress_sw_output_types'] = str_int_vals( 

955 self.stress_sw_output_types) 

956 d['str_tilt_sw_output_types'] = str_int_vals( 

957 self.tilt_sw_output_types) 

958 d['str_gravity_sw_output_types'] = str_int_vals( 

959 self.gravity_sw_output_types) 

960 

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

962#=============================================================================== 

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

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

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

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

967# read the accompanying READ.ME file. 

968# 

969# written by Rongjiang Wang 

970# GeoForschungsZentrum Potsdam 

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

972# phone +49 331 2881209 

973# fax +49 331 2881204 

974# 

975# Last modified: Potsdam, July, 2008 

976# 

977# References: 

978# 

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

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

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

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

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

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

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

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

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

988# International, 161, 191-196. 

989# 

990################################################################# 

991## ## 

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

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

994## ## 

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

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

997## downward. ## 

998## ## 

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

1000## ## 

1001################################################################# 

1002#=============================================================================== 

1003# OBSERVATION ARRAY 

1004# ================= 

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

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

1007# 

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

1009# 

1010# 2. number of positions: nrec 

1011# 

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

1013# 

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

1015# 

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

1017# 

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

1019# 

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

1021# 

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

1023# 

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

1025# 

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

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

1028# 

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

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

1031#=============================================================================== 

1032 %(sw_observation_type)i 

1033%(str_npoints)s 

1034%(str_observation)s 

1035#=============================================================================== 

1036# OUTPUTS 

1037# ======= 

1038# 

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

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

1041# 

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

1043# data: 

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

1045# 

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

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

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

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

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

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

1052# criterion) 

1053# 

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

1055# data: 

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

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

1058# oriented faults; 

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

1060# optimally oriented faults. 

1061# 

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

1063# 

1064# 3. output directory in char format: outdir 

1065# 

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

1067# 

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

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

1070# 

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

1072# 

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

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

1075# 

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

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

1078# 

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

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

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

1082# 

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

1084# as used for the Green's functions 

1085# 

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

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

1088# 

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

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

1091#13. ... 

1092# 

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

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

1095#=============================================================================== 

1096 %(sw_los_displacement)i %(str_los_vector)s 

1097 %(sw_coulomb_stress)i %(str_coulomb_master_field)s 

1098 %(str_pscmp_outdir)s 

1099 %(str_displ_sw_output_types)s 

1100 %(str_outdispl_filenames)s 

1101 %(str_stress_sw_output_types)s 

1102 %(str_outstress_filenames)s 

1103 %(str_tilt_sw_output_types)s %(str_gravity_sw_output_types)s 

1104 %(str_outtilt_filenames)s %(str_outgravity_filenames)s 

1105 %(n_snapshots)i 

1106%(str_times_snapshots)s 

1107#=============================================================================== 

1108# 

1109# GREEN'S FUNCTION DATABASE 

1110# ========================= 

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

1112# 

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

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

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

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

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

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

1119# 

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

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

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

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

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

1125# sources, respectively. 

1126# 

1127#=============================================================================== 

1128 %(str_psgrn_outdir)s 

1129 %(str_indispl_filenames)s 

1130 %(str_instress_filenames)s 

1131 %(str_intilt_filenames)s %(str_ingravity_filenames)s 

1132#=============================================================================== 

1133# RECTANGULAR SUBFAULTS 

1134# ===================== 

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

1136# 

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

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

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

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

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

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

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

1144# 

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

1146# 

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

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

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

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

1151# 

1152# 

1153# 3. ... for the 2. subfault ... 

1154# ... 

1155# N 

1156# / 

1157# /| strike 

1158# +------------------------ 

1159# |\ p . \ W 

1160# :-\ i . \ i 

1161# | \ l . \ d 

1162# :90 \ S . \ t 

1163# |-dip\ . \ h 

1164# : \. | rake \ 

1165# Z ------------------------- 

1166# L e n g t h 

1167# 

1168# Simulation of a Mogi source: 

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

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

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

1172# each other. 

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

1174# ratio at the source depth. 

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

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

1177#=============================================================================== 

1178# n_faults 

1179#------------------------------------------------------------------------------- 

1180 %(n_patches)i 

1181#------------------------------------------------------------------------------- 

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

1183# [-] [deg] [deg] [km] [km] [km] [deg] [deg] [-] [-] [day] 

1184# pos_s pos_d slp_stk slp_ddip open 

1185# [km] [km] [m] [m] [m] 

1186#------------------------------------------------------------------------------- 

1187%(patches_str)s 

1188#================================end of input=================================== 

1189''' # noqa 

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

1191 

1192 

1193class PsGrnPsCmpConfig(Object): 

1194 ''' 

1195 Combined config Object of PsGrn and PsCmp. 

1196 

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

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

1199 carefully. 

1200 ''' 

1201 tmin_days = Float.T( 

1202 default=0., 

1203 help='Min. time in days') 

1204 tmax_days = Float.T( 

1205 default=1., 

1206 help='Max. time in days') 

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

1208 

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

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

1211 

1212 

1213class PsCmpError(gf.store.StoreError): 

1214 pass 

1215 

1216 

1217class Interrupted(gf.store.StoreError): 

1218 def __str__(self): 

1219 return 'Interrupted.' 

1220 

1221 

1222class PsCmpRunner(object): 

1223 ''' 

1224 Wrapper object to execute the program fomosto_pscmp with the specified 

1225 configuration. 

1226 

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

1228 results are stored 

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

1230 ''' 

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

1232 if tmp is not None: 

1233 tmp = os.path.abspath(tmp) 

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

1235 self.keep_tmp = keep_tmp 

1236 self.config = None 

1237 

1238 def run(self, config): 

1239 ''' 

1240 Run the program! 

1241 

1242 :param config: :py:class:`PsCmpConfigFull` 

1243 ''' 

1244 self.config = config 

1245 

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

1247 

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

1249 input_str = config.string_for_config() 

1250 

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

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

1253 

1254 f.write(input_str) 

1255 

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

1257 

1258 old_wd = os.getcwd() 

1259 os.chdir(self.tempdir) 

1260 

1261 interrupted = [] 

1262 

1263 def signal_handler(signum, frame): 

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

1265 interrupted.append(True) 

1266 

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

1268 try: 

1269 try: 

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

1271 close_fds=True) 

1272 

1273 except OSError as err: 

1274 os.chdir(old_wd) 

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

1276 raise PsCmpError( 

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

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

1279on 

1280 

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

1282 

1283''' % program) 

1284 

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

1286 

1287 finally: 

1288 signal.signal(signal.SIGINT, original) 

1289 

1290 if interrupted: 

1291 raise KeyboardInterrupt() 

1292 

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

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

1295 

1296 errmsg = [] 

1297 if proc.returncode != 0: 

1298 errmsg.append( 

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

1300 

1301 if error_str: 

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

1303 

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

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

1306 

1307 if errmsg: 

1308 self.keep_tmp = True 

1309 

1310 os.chdir(old_wd) 

1311 raise PsCmpError(''' 

1312===== begin pscmp input ===== 

1313{pscmp_input}===== end pscmp input ===== 

1314===== begin pscmp output ===== 

1315{pscmp_output}===== end pscmp output ===== 

1316===== begin pscmp error ===== 

1317{pscmp_error}===== end pscmp error ===== 

1318{error_messages} 

1319pscmp has been invoked as "{call}" 

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

1321 pscmp_input=input_str, 

1322 pscmp_output=output_str, 

1323 pscmp_error=error_str, 

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

1325 call=program, 

1326 dir=self.tempdir) 

1327 .lstrip()) 

1328 

1329 self.pscmp_output = output_str 

1330 self.pscmp_error = error_str 

1331 

1332 os.chdir(old_wd) 

1333 

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

1335 ''' 

1336 Get the resulting components from the stored directory. 

1337 Be careful: The z-component is downward positive! 

1338 

1339 :param component: string, the component to retrieve from the 

1340 result directory, may be: 

1341 "displ": displacement, n x 3 array 

1342 "stress": stresses n x 6 array 

1343 "tilt': tilts n x 3 array, 

1344 "gravity': gravity n x 2 array 

1345 "all": all the above together 

1346 ''' 

1347 assert self.config.snapshots is not None 

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

1349 

1350 output = [] 

1351 for fn in fns: 

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

1353 continue 

1354 

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

1356 

1357 try: 

1358 _, idxs = pscmp_component_mapping[component] 

1359 except KeyError: 

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

1361 component, ', '.join( 

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

1363 

1364 istart, iend = idxs 

1365 

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

1367 

1368 return output 

1369 

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

1371 ''' 

1372 Load snapshot data array and return specified components. 

1373 Transform array component and receiver wise to list of 

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

1375 ''' 

1376 

1377 distances = self.config.observation.distances 

1378 

1379 output_list = self.get_results(component=component) 

1380 deltat = self.config.snapshots.deltat 

1381 

1382 nrec, ncomp = output_list[0].shape 

1383 

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

1385 

1386 comps, _ = pscmp_component_mapping[component] 

1387 

1388 outtraces = [] 

1389 for row in range(nrec): 

1390 for col, comp in enumerate(comps): 

1391 tr = trace.Trace( 

1392 '', '%04i' % row, '', comp, 

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

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

1395 outtraces.append(tr) 

1396 

1397 return outtraces 

1398 

1399 def __del__(self): 

1400 if self.tempdir: 

1401 if not self.keep_tmp: 

1402 shutil.rmtree(self.tempdir) 

1403 self.tempdir = None 

1404 else: 

1405 logger.warning( 

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

1407 

1408 

1409class PsGrnCmpGFBuilder(gf.builder.Builder): 

1410 nsteps = 2 

1411 

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

1413 force=False): 

1414 

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

1416 

1417 storeconf = self.store.config 

1418 

1419 dummy_lat = 10.0 

1420 dummy_lon = 10.0 

1421 

1422 depths = storeconf.coords[0] 

1423 lats = num.ones_like(depths) * dummy_lat 

1424 lons = num.ones_like(depths) * dummy_lon 

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

1426 

1427 self.shear_moduli = storeconf.get_shear_moduli( 

1428 lat=dummy_lat, 

1429 lon=dummy_lon, 

1430 points=points, 

1431 interpolation='multilinear') 

1432 

1433 self.lambda_moduli = storeconf.get_lambda_moduli( 

1434 lat=dummy_lat, 

1435 lon=dummy_lon, 

1436 points=points, 

1437 interpolation='multilinear') 

1438 

1439 if step == 0: 

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

1441 else: 

1442 if block_size is None: 

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

1444 

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

1446 block_size = block_size[1:] 

1447 

1448 gf.builder.Builder.__init__( 

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

1450 

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

1452 

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

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

1455 deltat_days = 1. / storeconf.sample_rate / day 

1456 

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

1458 cg.max_time = baseconf.tmax_days 

1459 else: 

1460 warnings.warn( 

1461 'PsGrnConfig is defining n_times_samples and max_time,' 

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

1463 FutureWarning) 

1464 

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

1466 if cc.snapshots is None: 

1467 deltat_days = 1. / storeconf.sample_rate / day 

1468 

1469 cc.snapshots = PsCmpSnapshots( 

1470 tmin=baseconf.tmin_days, 

1471 tmax=baseconf.tmax_days, 

1472 deltatdays=deltat_days) 

1473 else: 

1474 warnings.warn( 

1475 'PsCmpConfig is defining snapshots,' 

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

1477 FutureWarning) 

1478 

1479 cg.earthmodel_1d = storeconf.earthmodel_1d 

1480 

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

1482 

1483 cg.psgrn_outdir = gf_outpath + '/' 

1484 cc.psgrn_outdir = gf_outpath + '/' 

1485 

1486 util.ensuredir(gf_outpath) 

1487 

1488 self.cg = cg 

1489 self.cc = cc 

1490 

1491 def cleanup(self): 

1492 self.store.close() 

1493 

1494 def work_block(self, iblock): 

1495 

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

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

1498 self.get_block_extents(iblock) 

1499 mu = self.shear_moduli[iblock] 

1500 lame_lambda = self.lambda_moduli[iblock] 

1501 

1502 rz = self.store.config.receiver_depth 

1503 else: 

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

1505 self.get_block_extents(iblock) 

1506 

1507 fc = self.store.config 

1508 cc = copy.deepcopy(self.cc) 

1509 cg = copy.deepcopy(self.cg) 

1510 

1511 dx = fc.distance_delta 

1512 

1513 logger.info( 

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

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

1516 

1517 if self.step == 0: 

1518 if cg.gf_depth_spacing is None: 

1519 gf_depth_spacing = fc.source_depth_delta 

1520 else: 

1521 gf_depth_spacing = cg.gf_depth_spacing * km 

1522 

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

1524 gf_depth_spacing) + 1 

1525 

1526 if cg.gf_distance_spacing is None: 

1527 gf_distance_spacing = fc.distance_delta 

1528 else: 

1529 gf_distance_spacing = cg.gf_distance_spacing * km 

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

1531 gf_distance_spacing) + 1 

1532 

1533 cg.depth_grid = PsGrnSpatialSampling( 

1534 n_steps=n_steps_depth, 

1535 start_distance=fc.source_depth_min / km, 

1536 end_distance=fc.source_depth_max / km) 

1537 

1538 cg.distance_grid = PsGrnSpatialSampling( 

1539 n_steps=n_steps_distance, 

1540 start_distance=fc.distance_min / km, 

1541 end_distance=fc.distance_max / km) 

1542 

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

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

1545 

1546 else: 

1547 distances = num.linspace( 

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

1549 

1550 # fomosto sample rate in s, pscmp takes days 

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

1552 cc.snapshots = PsCmpSnapshots( 

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

1554 cc.observation = PsCmpProfile( 

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

1556 slon=0.0, 

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

1558 elon=0.0, 

1559 n_steps=len(distances), 

1560 distances=distances) 

1561 

1562 runner = PsCmpRunner(keep_tmp=False) 

1563 

1564 mtsize = float(dx * cc.rectangular_fault_size_factor) 

1565 

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

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

1568 

1569 mui = 1. / mu 

1570 

1571 gfmapping = [ 

1572 (('nn',), 

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

1574 (('ne',), 

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

1576 (('nd',), 

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

1578 (('ed',), 

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

1580 (('dd',), 

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

1582 (('ee',), 

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

1584 ] 

1585 

1586 for mt, gfmap in gfmapping: 

1587 cc.rectangular_source_patches = [] 

1588 for idx in mt: 

1589 pmt = PsCmpMomentTensor( 

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

1591 lon=0.0, 

1592 depth=float(sz), 

1593 width=mtsize, 

1594 length=mtsize, 

1595 idx=idx) 

1596 

1597 cc.rectangular_source_patches.extend( 

1598 pmt.to_rfs(nullf)) 

1599 

1600 runner.run(cc) 

1601 

1602 rawtraces = runner.get_traces() 

1603 

1604 interrupted = [] 

1605 

1606 def signal_handler(signum, frame): 

1607 interrupted.append(True) 

1608 

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

1610 self.store.lock() 

1611 duplicate_inserts = 0 

1612 

1613 try: 

1614 for itr, tr in enumerate(rawtraces): 

1615 if tr.channel in gfmap: 

1616 x = tr.meta['distance'] 

1617 ig, factor = gfmap[tr.channel] 

1618 

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

1620 args = (sz, x, ig) 

1621 else: 

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

1623 

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

1625 

1626 gf_tr.data *= factor 

1627 

1628 try: 

1629 self.store.put(args, gf_tr) 

1630 except gf.store.DuplicateInsert: 

1631 duplicate_inserts += 1 

1632 

1633 finally: 

1634 if duplicate_inserts: 

1635 logger.warning( 

1636 '%i insertions skipped (duplicates)' 

1637 % duplicate_inserts) 

1638 

1639 self.store.unlock() 

1640 signal.signal(signal.SIGINT, original) 

1641 

1642 if interrupted: 

1643 raise KeyboardInterrupt() 

1644 

1645 logger.info( 

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

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

1648 

1649 

1650def init(store_dir, variant): 

1651 if variant is None: 

1652 variant = '2008a' 

1653 

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

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

1656 

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

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

1659 

1660 c = PsGrnPsCmpConfig() 

1661 

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

1663 

1664 # Initialising a viscous mantle 

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

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

1667 mantle.burger_eta1 = 5e17 

1668 mantle.burger_eta2 = 1e19 

1669 mantle.burger_alpha = 1. 

1670 

1671 config = gf.meta.ConfigTypeA( 

1672 id=store_id, 

1673 ncomponents=10, 

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

1675 receiver_depth=0. * km, 

1676 source_depth_min=0. * km, 

1677 source_depth_max=15. * km, 

1678 source_depth_delta=.5 * km, 

1679 distance_min=0. * km, 

1680 distance_max=50. * km, 

1681 distance_delta=1. * km, 

1682 earthmodel_1d=cake_mod, 

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

1684 tabulated_phases=[]) # dummy list 

1685 

1686 c.validate() 

1687 config.validate() 

1688 return gf.store.Store.create_editables( 

1689 store_dir, 

1690 config=config, 

1691 extra={'psgrn_pscmp': c}) 

1692 

1693 

1694def build(store_dir, 

1695 force=False, 

1696 nworkers=None, 

1697 continue_=False, 

1698 step=None, 

1699 iblock=None): 

1700 

1701 return PsGrnCmpGFBuilder.build( 

1702 store_dir, force=force, nworkers=nworkers, continue_=continue_, 

1703 step=step, iblock=iblock)