1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import logging
8import warnings
9import os
10from os.path import join as pjoin
11import signal
12import shutil
13import copy
15import math
16import numpy as num
18from tempfile import mkdtemp
19from subprocess import Popen, PIPE
21from pyrocko.guts import Float, Int, Tuple, List, Object, String
22from pyrocko.model import Location
23from pyrocko import gf, util, trace, cake
26# how to call the programs
27program_bins = {
28 'pscmp.2008a': 'fomosto_pscmp2008a',
29 'psgrn.2008a': 'fomosto_psgrn2008a'
30}
32psgrn_displ_names = ('uz', 'ur', 'ut')
33psgrn_stress_names = ('szz', 'srr', 'stt', 'szr', 'srt', 'stz')
34psgrn_tilt_names = ('tr', 'tt', 'rot')
35psgrn_gravity_names = ('gd', 'gr')
36psgrn_components = 'ep ss ds cl'.split()
38km = 1000.
39day = 3600. * 24
41guts_prefix = 'pf'
42logger = logging.getLogger('pyrocko.fomosto.psgrn_pscmp')
45def have_backend():
46 for cmd in [[exe] for exe in program_bins.values()]:
47 try:
48 p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE)
49 (stdout, stderr) = p.communicate()
51 except OSError:
52 return False
54 return True
57def nextpow2(i):
58 return 2 ** int(math.ceil(math.log(i) / math.log(2.)))
61def str_float_vals(vals):
62 return ' '.join('%e' % val for val in vals)
65def str_int_vals(vals):
66 return ' '.join('%i' % val for val in vals)
69def str_str_vals(vals):
70 return ' '.join("'%s'" % val for val in vals)
73def cake_model_to_config(mod):
74 srows = []
75 for ir, row in enumerate(mod.to_scanlines(get_burgers=True)):
76 depth, vp, vs, rho, qp, qs, eta1, eta2, alpha = row
77 # replace qs with etas = 0.
78 row = [depth / km, vp / km, vs / km, rho, eta1, eta2, alpha]
79 srows.append('%i %15s' % (ir + 1, str_float_vals(row)))
81 return '\n'.join(srows), len(srows)
84class PsGrnSpatialSampling(Object):
86 '''
87 Definition of spatial sampling for PSGRN.
89 Note: attributes in this class use non-standard units [km] to be consistent
90 with PSGRN text file input. Please read the documentation carefully.
91 '''
93 n_steps = Int.T(default=10)
94 start_distance = Float.T(default=0.) # start sampling [km] from source
95 end_distance = Float.T(default=100.) # end
97 def string_for_config(self):
98 return '%i %15e %15e' % (self.n_steps, self.start_distance,
99 self.end_distance)
102class PsGrnConfig(Object):
104 '''
105 Configuration for PSGRN.
107 Note: attributes in this class use non-standard units [km] and [days] to
108 be consistent with PSGRN text file input. Please read the documentation
109 carefully.
110 '''
112 version = String.T(default='2008a')
114 sampling_interval = Float.T(
115 default=1.0,
116 help='Exponential sampling 1.- equidistant, > 1. decreasing sampling'
117 ' with distance')
118 n_time_samples = Int.T(
119 optional=True,
120 help='Number of temporal GF samples up to max_time. Has to be equal'
121 ' to a power of 2! If not, next power of 2 is taken.')
122 max_time = Float.T(
123 optional=True,
124 help='Maximum time [days] after seismic event.')
126 gf_depth_spacing = Float.T(
127 optional=True,
128 help='Depth spacing [km] for the observation points. '
129 'If not defined depth spacing from the `StoreConfig` is used')
130 gf_distance_spacing = Float.T(
131 optional=True,
132 help='Spatial spacing [km] for the observation points. '
133 'If not defined distance spacing from the `StoreConfig` is used')
134 observation_depth = Float.T(
135 default=0.,
136 help='Depth of observation points [km]')
138 def items(self):
139 return dict(self.T.inamevals(self))
142class PsGrnConfigFull(PsGrnConfig):
144 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True)
145 psgrn_outdir = String.T(default='psgrn_green/')
147 distance_grid = PsGrnSpatialSampling.T(PsGrnSpatialSampling.D())
148 depth_grid = PsGrnSpatialSampling.T(PsGrnSpatialSampling.D())
150 sw_source_regime = Int.T(default=1) # 1-continental, 0-ocean
151 sw_gravity = Int.T(default=0)
153 accuracy_wavenumber_integration = Float.T(default=0.025)
155 displ_filenames = Tuple.T(3, String.T(), default=psgrn_displ_names)
156 stress_filenames = Tuple.T(6, String.T(), default=psgrn_stress_names)
157 tilt_filenames = Tuple.T(3, String.T(), psgrn_tilt_names)
158 gravity_filenames = Tuple.T(2, String.T(), psgrn_gravity_names)
160 @staticmethod
161 def example():
162 conf = PsGrnConfigFull()
163 conf.earthmodel_1d = cake.load_model().extract(depth_max=100*km)
164 conf.psgrn_outdir = 'TEST_psgrn_functions/'
165 return conf
167 def string_for_config(self):
169 assert self.earthmodel_1d is not None
171 d = self.__dict__.copy()
173 model_str, nlines = cake_model_to_config(self.earthmodel_1d)
174 d['n_t2'] = nextpow2(self.n_time_samples)
175 d['n_model_lines'] = nlines
176 d['model_lines'] = model_str
178 d['str_psgrn_outdir'] = "'%s'" % './'
180 d['str_displ_filenames'] = str_str_vals(self.displ_filenames)
181 d['str_stress_filenames'] = str_str_vals(self.stress_filenames)
182 d['str_tilt_filenames'] = str_str_vals(self.tilt_filenames)
183 d['str_gravity_filenames'] = str_str_vals(self.gravity_filenames)
185 d['str_distance_grid'] = self.distance_grid.string_for_config()
186 d['str_depth_grid'] = self.depth_grid.string_for_config()
188 template = '''# autogenerated PSGRN input by psgrn.py
189#=============================================================================
190# This is input file of FORTRAN77 program "psgrn08a" for computing responses
191# (Green's functions) of a multi-layered viscoelastic halfspace to point
192# dislocation sources buried at different depths. All results will be stored in
193# the given directory and provide the necessary data base for the program
194# "pscmp07a" for computing time-dependent deformation, geoid and gravity changes
195# induced by an earthquake with extended fault planes via linear superposition.
196# For more details, please read the accompanying READ.ME file.
197#
198# written by Rongjiang Wang
199# GeoForschungsZentrum Potsdam
200# e-mail: wang@gfz-potsdam.de
201# phone +49 331 2881209
202# fax +49 331 2881204
203#
204# Last modified: Potsdam, Jan, 2008
205#
206#################################################################
207## ##
208## Cylindrical coordinates (Z positive downwards!) are used. ##
209## ##
210## If not specified otherwise, SI Unit System is used overall! ##
211## ##
212#################################################################
213#
214#------------------------------------------------------------------------------
215#
216# PARAMETERS FOR SOURCE-OBSERVATION CONFIGURATIONS
217# ================================================
218# 1. the uniform depth of the observation points [km], switch for oceanic (0)
219# or continental(1) earthquakes;
220# 2. number of (horizontal) observation distances (> 1 and <= nrmax defined in
221# psgglob.h), start and end distances [km], ratio (>= 1.0) between max. and
222# min. sampling interval (1.0 for equidistant sampling);
223# 3. number of equidistant source depths (>= 1 and <= nzsmax defined in
224# psgglob.h), start and end source depths [km];
225#
226# r1,r2 = minimum and maximum horizontal source-observation
227# distances (r2 > r1).
228# zs1,zs2 = minimum and maximum source depths (zs2 >= zs1 > 0).
229#
230# Note that the same sampling rates dr_min and dzs will be used later by the
231# program "pscmp07a" for discretizing the finite source planes to a 2D grid
232# of point sources.
233#------------------------------------------------------------------------------
234 %(observation_depth)e %(sw_source_regime)i
235 %(str_distance_grid)s %(sampling_interval)e
236 %(str_depth_grid)s
237#------------------------------------------------------------------------------
238#
239# PARAMETERS FOR TIME SAMPLING
240# ============================
241# 1. number of time samples (<= ntmax def. in psgglob.h) and time window [days].
242#
243# Note that nt (> 0) should be power of 2 (the fft-rule). If nt = 1, the
244# coseismic (t = 0) changes will be computed; If nt = 2, the coseismic
245# (t = 0) and steady-state (t -> infinity) changes will be computed;
246# Otherwise, time series for the given time samples will be computed.
247#
248#------------------------------------------------------------------------------
249 %(n_t2)i %(max_time)f
250#------------------------------------------------------------------------------
251#
252# PARAMETERS FOR WAVENUMBER INTEGRATION
253# =====================================
254# 1. relative accuracy of the wave-number integration (suggested: 0.1 - 0.01)
255# 2. factor (> 0 and < 1) for including influence of earth's gravity on the
256# deformation field (e.g. 0/1 = without / with 100percent gravity effect).
257#------------------------------------------------------------------------------
258 %(accuracy_wavenumber_integration)e
259 %(sw_gravity)i
260#------------------------------------------------------------------------------
261#
262# PARAMETERS FOR OUTPUT FILES
263# ===========================
264#
265# 1. output directory
266# 2. file names for 3 displacement components (uz, ur, ut)
267# 3. file names for 6 stress components (szz, srr, stt, szr, srt, stz)
268# 4. file names for radial and tangential tilt components (as measured by a
269# borehole tiltmeter), rigid rotation of horizontal plane, geoid and gravity
270# changes (tr, tt, rot, gd, gr)
271#
272# Note that all file or directory names should not be longer than 80
273# characters. Directory and subdirectoy names must be separated and ended
274# by / (unix) or \ (dos)! All file names should be given without extensions
275# that will be appended automatically by ".ep" for the explosion (inflation)
276# source, ".ss" for the strike-slip source, ".ds" for the dip-slip source,
277# and ".cl" for the compensated linear vector dipole source)
278#
279#------------------------------------------------------------------------------
280 %(str_psgrn_outdir)s
281 %(str_displ_filenames)s
282 %(str_stress_filenames)s
283 %(str_tilt_filenames)s %(str_gravity_filenames)s
284#------------------------------------------------------------------------------
285#
286# GLOBAL MODEL PARAMETERS
287# =======================
288# 1. number of data lines of the layered model (<= lmax as defined in psgglob.h)
289#
290# The surface and the upper boundary of the half-space as well as the
291# interfaces at which the viscoelastic parameters are continuous, are all
292# defined by a single data line; All other interfaces, at which the
293# viscoelastic parameters are discontinuous, are all defined by two data
294# lines (upper-side and lower-side values). This input format could also be
295# used for a graphic plot of the layered model. Layers which have different
296# parameter values at top and bottom, will be treated as layers with a
297# constant gradient, and will be discretised to a number of homogeneous
298# sublayers. Errors due to the discretisation are limited within about
299# 5percent (changeable, see psgglob.h).
300#
301# 2.... parameters of the multilayered model
302#
303# Burgers rheology (a Kelvin-Voigt body and a Maxwell body in series
304# connection) for relaxation of shear modulus is implemented. No relaxation
305# of compressional modulus is considered.
306#
307# eta1 = transient viscosity (dashpot of the Kelvin-Voigt body; <= 0 means
308# infinity value)
309# eta2 = steady-state viscosity (dashpot of the Maxwell body; <= 0 means
310# infinity value)
311# alpha = ratio between the effective and the unrelaxed shear modulus
312# = mu1/(mu1+mu2) (> 0 and <= 1)
313#
314# Special cases:
315# (1) Elastic: eta1 and eta2 <= 0 (i.e. infinity); alpha meaningless
316# (2) Maxwell body: eta1 <= 0 (i.e. eta1 = infinity)
317# or alpha = 1 (i.e. mu1 = infinity)
318# (3) Standard-Linear-Solid: eta2 <= 0 (i.e. infinity)
319#------------------------------------------------------------------------------
320 %(n_model_lines)i |int: no_model_lines;
321#------------------------------------------------------------------------------
322# no depth[km] vp[km/s] vs[km/s] rho[kg/m^3] eta1[Pa*s] eta2[Pa*s] alpha
323#------------------------------------------------------------------------------
324%(model_lines)s
325#=======================end of input===========================================
326''' # noqa
327 return (template % d).encode('ascii')
330class PsGrnError(gf.store.StoreError):
331 pass
334def remove_if_exists(fn, force=False):
335 if os.path.exists(fn):
336 if force:
337 os.remove(fn)
338 else:
339 raise gf.CannotCreate('file %s already exists' % fn)
342class PsGrnRunner(object):
343 '''
344 Wrapper object to execute the program fomosto_psgrn.
345 '''
347 def __init__(self, outdir):
348 outdir = os.path.abspath(outdir)
349 if not os.path.exists(outdir):
350 os.mkdir(outdir)
351 self.outdir = outdir
352 self.config = None
354 def run(self, config, force=False):
355 '''
356 Run the program with the specified configuration.
358 :param config: :py:class:`PsGrnConfigFull`
359 :param force: boolean, set true to overwrite existing output
360 '''
361 self.config = config
363 input_fn = pjoin(self.outdir, 'input')
365 remove_if_exists(input_fn, force=force)
367 with open(input_fn, 'wb') as f:
368 input_str = config.string_for_config()
369 logger.debug('===== begin psgrn input =====\n'
370 '%s===== end psgrn input =====' % input_str.decode())
371 f.write(input_str)
372 program = program_bins['psgrn.%s' % config.version]
374 old_wd = os.getcwd()
376 os.chdir(self.outdir)
378 interrupted = []
380 def signal_handler(signum, frame):
381 os.kill(proc.pid, signal.SIGTERM)
382 interrupted.append(True)
384 original = signal.signal(signal.SIGINT, signal_handler)
385 try:
386 try:
387 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
389 except OSError:
390 os.chdir(old_wd)
391 raise PsGrnError(
392 '''could not start psgrn executable: "%s"
393Available fomosto backends and download links to the modelling codes are listed
394on
396 https://pyrocko.org/docs/current/apps/fomosto/backends.html
398''' % program)
400 (output_str, error_str) = proc.communicate(b'input\n')
402 finally:
403 signal.signal(signal.SIGINT, original)
405 if interrupted:
406 raise KeyboardInterrupt()
408 logger.debug('===== begin psgrn output =====\n'
409 '%s===== end psgrn output =====' % output_str.decode())
411 errmess = []
412 if proc.returncode != 0:
413 errmess.append(
414 'psgrn had a non-zero exit state: %i' % proc.returncode)
416 if error_str:
417 logger.warning(
418 'psgrn emitted something via stderr: \n\n%s'
419 % error_str.decode())
420 # errmess.append('psgrn emitted something via stderr')
422 if output_str.lower().find(b'error') != -1:
423 errmess.append("the string 'error' appeared in psgrn output")
425 if errmess:
426 os.chdir(old_wd)
427 raise PsGrnError('''
428===== begin psgrn input =====
429%s===== end psgrn input =====
430===== begin psgrn output =====
431%s===== end psgrn output =====
432===== begin psgrn error =====
433%s===== end psgrn error =====
434%s
435psgrn has been invoked as "%s"
436in the directory %s'''.lstrip() % (
437 input_str.decode(), output_str.decode(), error_str.decode(),
438 '\n'.join(errmess), program, self.outdir))
440 self.psgrn_output = output_str
441 self.psgrn_error = error_str
443 os.chdir(old_wd)
446pscmp_displ_names = ('un', 'ue', 'ud')
447pscmp_stress_names = ('snn', 'see', 'sdd', 'sne', 'snd', 'sed')
448pscmp_tilt_names = ('tn', 'te', 'rot')
449pscmp_gravity_names = ('gd', 'gr')
450pscmp_all_names = pscmp_displ_names + pscmp_stress_names + pscmp_tilt_names + \
451 pscmp_gravity_names
453pscmp_component_mapping = {
454 'displ': (pscmp_displ_names, (2, 5)),
455 'stress': (pscmp_stress_names, (5, 11)),
456 'tilt': (pscmp_tilt_names, (11, 14)),
457 'gravity': (pscmp_gravity_names, (14, 16)),
458 'all': (pscmp_all_names, (2, 16)),
459 }
462def dsin(value):
463 return num.sin(value * num.pi / 180.)
466def dcos(value):
467 return num.cos(value * num.pi / 180.)
470def distributed_fault_patches_to_config(patches):
471 '''
472 Input: List of PsCmpRectangularSource(s)
473 '''
474 srows = []
475 for i, patch in enumerate(patches):
476 srows.append('%i %s' % (i + 1, patch.string_for_config()))
478 return '\n'.join(srows), len(srows)
481class PsCmpObservation(Object):
482 pass
485class PsCmpScatter(PsCmpObservation):
486 '''
487 Scattered observation points.
488 '''
489 lats = List.T(Float.T(), optional=True, default=[10.4, 10.5])
490 lons = List.T(Float.T(), optional=True, default=[12.3, 13.4])
492 def string_for_config(self):
493 srows = []
494 for lat, lon in zip(self.lats, self.lons):
495 srows.append('(%15f, %15f)' % (lat, lon))
497 self.sw = 0
498 return ' %i' % (len(srows)), '\n'.join(srows)
501class PsCmpProfile(PsCmpObservation):
502 '''
503 Calculation along observation profile.
504 '''
505 n_steps = Int.T(default=10)
506 slat = Float.T(
507 default=0.,
508 help='Profile start latitude')
509 slon = Float.T(
510 default=0.,
511 help='Profile start longitude')
512 elat = Float.T(
513 default=0.,
514 help='Profile end latitude')
515 elon = Float.T(
516 default=0.,
517 help='Profile end longitude')
518 distances = List.T(
519 Float.T(),
520 optional=True,
521 help='Distances [m] for each point on profile from start to end.')
523 def string_for_config(self):
524 self.sw = 1
526 return ' %i' % (self.n_steps), \
527 ' ( %15f, %15f ), ( %15f, %15f )' % (
528 self.slat, self.slon, self.elat, self.elon)
531class PsCmpArray(PsCmpObservation):
532 '''
533 Calculation on a grid.
534 '''
535 n_steps_lat = Int.T(default=10)
536 n_steps_lon = Int.T(default=10)
537 slat = Float.T(
538 default=0.,
539 help='Array start latitude')
540 slon = Float.T(
541 default=0.,
542 help='Array start longitude')
543 elat = Float.T(
544 default=0.,
545 help='Array end latitude')
546 elon = Float.T(
547 default=0.,
548 help='Array end longitude')
550 def string_for_config(self):
551 self.sw = 2
553 return ' %i %15f %15f ' % (
554 self.n_steps_lat, self.slat, self.elat), \
555 ' %i %15f %15f ' % (
556 self.n_steps_lon, self.slon, self.elon)
559class PsCmpRectangularSource(Location, gf.seismosizer.Cloneable):
560 '''
561 Rectangular Source for the input geometry of the active fault.
563 Input parameters have to be in:
564 [deg] for reference point (lat, lon) and angles (rake, strike, dip)
565 [m] shifting with respect to reference position
566 [m] for fault dimensions and source depth. The default shift of the
567 origin (:py:attr`pos_s`, :py:attr:`pos_d`) with respect to the reference
568 coordinates
569 (lat, lon) is zero, which implies that the reference is the center of
570 the fault plane!
571 The calculation point is always the center of the fault-plane!
572 Setting :py:attr`pos_s` or :py:attr`pos_d` moves the fault point with
573 respect to the origin along strike and dip direction, respectively!
574 '''
575 length = Float.T(default=6.0 * km)
576 width = Float.T(default=5.0 * km)
577 strike = Float.T(default=0.0)
578 dip = Float.T(default=90.0)
579 rake = Float.T(default=0.0)
580 torigin = Float.T(default=0.0)
582 slip = Float.T(optional=True, default=1.0)
584 pos_s = Float.T(optional=True, default=None)
585 pos_d = Float.T(optional=True, default=None)
586 opening = Float.T(default=0.0)
588 def update(self, **kwargs):
589 '''
590 Change some of the source models parameters.
592 Example::
594 >>> from pyrocko import gf
595 >>> s = gf.DCSource()
596 >>> s.update(strike=66., dip=33.)
597 >>> print(s)
598 --- !pf.DCSource
599 depth: 0.0
600 time: 1970-01-01 00:00:00
601 magnitude: 6.0
602 strike: 66.0
603 dip: 33.0
604 rake: 0.0
606 '''
607 for (k, v) in kwargs.items():
608 self[k] = v
610 @property
611 def dip_slip(self):
612 return float(self.slip * dsin(self.rake) * (-1))
614 @property
615 def strike_slip(self):
616 return float(self.slip * dcos(self.rake))
618 def string_for_config(self):
620 if self.pos_s or self.pos_d is None:
621 self.pos_s = 0.
622 self.pos_d = 0.
624 tempd = copy.deepcopy(self.__dict__)
625 tempd['dip_slip'] = self.dip_slip
626 tempd['strike_slip'] = self.strike_slip
627 tempd['effective_lat'] = self.effective_lat
628 tempd['effective_lon'] = self.effective_lon
629 tempd['depth'] /= km
630 tempd['length'] /= km
631 tempd['width'] /= km
633 return '%(effective_lat)15f %(effective_lon)15f %(depth)15f' \
634 '%(length)15f %(width)15f %(strike)15f' \
635 '%(dip)15f 1 1 %(torigin)15f \n %(pos_s)15f %(pos_d)15f ' \
636 '%(strike_slip)15f %(dip_slip)15f %(opening)15f' % tempd
639MTIso = {
640 'nn': dict(strike=90., dip=90., rake=0., slip=0., opening=1.),
641 'ee': dict(strike=0., dip=90., rake=0., slip=0., opening=1.),
642 'dd': dict(strike=0., dip=0., rake=-90., slip=0., opening=1.),
643 }
645MTDev = {
646 'ne': dict(strike=90., dip=90., rake=180., slip=1., opening=0.),
647 'nd': dict(strike=180., dip=0., rake=0., slip=1., opening=0.),
648 'ed': dict(strike=270., dip=0., rake=0., slip=1., opening=0.),
649 }
652def get_nullification_factor(mu, lame_lambda):
653 '''
654 Factor that needs to be multiplied to 2 of the tensile sources to
655 eliminate two of the isotropic components.
656 '''
657 return -lame_lambda / (2. * mu + 2. * lame_lambda)
660def get_trace_normalization_factor(mu, lame_lambda, nullification):
661 '''
662 Factor to be multiplied to elementary GF trace to have unit displacement.
663 '''
664 return 1. / ((2. * mu) + lame_lambda + (2. * lame_lambda * nullification))
667def get_iso_scaling_factors(mu, lame_lambda):
668 nullf = get_nullification_factor(mu, lame_lambda)
669 scale = get_trace_normalization_factor(mu, lame_lambda, nullf)
670 return nullf, scale
673class PsCmpTensileSF(Location, gf.seismosizer.Cloneable):
674 '''
675 Compound dislocation of 3 perpendicular, rectangular sources to approximate
676 an opening single force couple. NED coordinate system!
678 Warning: Mixed standard [m] / non-standard [days] units are used in this
679 class. Please read the documentation carefully.
680 '''
682 length = Float.T(
683 default=1.0 * km,
684 help='Length of source rectangle [m].')
685 width = Float.T(
686 default=1.0 * km,
687 help='Width of source rectangle [m].')
688 torigin = Float.T(
689 default=0.0,
690 help='Origin time [days]')
691 idx = String.T(
692 default='nn',
693 help='Axis index for opening direction; "nn", "ee", or "dd"')
695 def to_rfs(self, nullification):
697 cmpd = []
698 for comp, mt in MTIso.items():
699 params = copy.deepcopy(mt)
701 if comp != self.idx:
702 params = copy.deepcopy(mt)
703 params['opening'] *= nullification
705 kwargs = copy.deepcopy(self.__dict__)
706 kwargs.update(params)
707 kwargs.pop('idx')
708 kwargs.pop('_latlon')
709 cmpd.append(PsCmpRectangularSource(**kwargs))
711 return cmpd
714class PsCmpShearSF(Location, gf.seismosizer.Cloneable):
715 '''
716 Shear fault source model component.
718 Warning: Mixed standard [m] / non-standard [days] units are used in this
719 class. Please read the documentation carefully.
720 '''
722 length = Float.T(
723 default=1.0 * km,
724 help='Length of source rectangle [m].')
725 width = Float.T(
726 default=1.0 * km,
727 help='Width of source rectangle [m].')
728 torigin = Float.T(
729 default=0.0,
730 help='Origin time [days]')
731 idx = String.T(
732 default='ne',
733 help='Axis index for opening direction; "ne", "nd", or "ed"')
735 def to_rfs(self):
736 kwargs = copy.deepcopy(self.__dict__)
737 kwargs.pop('idx')
738 kwargs.pop('_latlon')
739 kwargs.update(MTDev[self.idx])
740 return [PsCmpRectangularSource(**kwargs)]
743class PsCmpMomentTensor(Location, gf.seismosizer.Cloneable):
744 '''
745 Mapping of Moment Tensor components to rectangular faults.
746 Only one component at a time valid! NED coordinate system!
748 Warning: Mixed standard [m] / non-standard [days] units are used in this
749 class. Please read the documentation carefully.
750 '''
751 length = Float.T(
752 default=1.0 * km,
753 help='Length of source rectangle [m].')
754 width = Float.T(
755 default=1.0 * km,
756 help='Width of source rectangle [m].')
757 torigin = Float.T(
758 default=0.0,
759 help='Origin time [days]')
760 idx = String.T(
761 default='nn',
762 help='Axis index for MT component; '
763 '"nn", "ee", "dd", "ne", "nd", or "ed".')
765 def to_rfs(self, nullification=-0.25):
766 kwargs = copy.deepcopy(self.__dict__)
767 kwargs.pop('_latlon')
769 if self.idx in MTIso:
770 tsf = PsCmpTensileSF(**kwargs)
771 return tsf.to_rfs(nullification)
773 elif self.idx in MTDev:
774 ssf = PsCmpShearSF(**kwargs)
775 return ssf.to_rfs()
777 else:
778 raise Exception('MT component not supported!')
781class PsCmpCoulombStress(Object):
782 pass
785class PsCmpCoulombStressMasterFault(PsCmpCoulombStress):
786 friction = Float.T(default=0.7)
787 skempton_ratio = Float.T(default=0.0)
788 master_fault_strike = Float.T(default=300.)
789 master_fault_dip = Float.T(default=15.)
790 master_fault_rake = Float.T(default=90.)
791 sigma1 = Float.T(default=1.e6, help='[Pa]')
792 sigma2 = Float.T(default=-1.e6, help='[Pa]')
793 sigma3 = Float.T(default=0.0, help='[Pa]')
795 def string_for_config(self):
796 return '%(friction)15e %(skempton_ratio)15e %(master_fault_strike)15f'\
797 '%(master_fault_dip)15f %(master_fault_rake)15f'\
798 '%(sigma1)15e %(sigma2)15e %(sigma3)15e' % self.__dict__
801class PsCmpSnapshots(Object):
802 '''
803 Snapshot time series definition.
805 Note: attributes in this class use non-standard units [days] to be
806 consistent with PSCMP text file input. Please read the documentation
807 carefully.
808 '''
810 tmin = Float.T(
811 default=0.0,
812 help='Time [days] after source time to start temporal sample'
813 ' snapshots.')
814 tmax = Float.T(
815 default=1.0,
816 help='Time [days] after source time to end temporal sample f.')
817 deltatdays = Float.T(
818 default=1.0,
819 help='Sample period [days].')
821 @property
822 def times(self):
823 nt = int(num.ceil((self.tmax - self.tmin) / self.deltatdays))
824 return num.linspace(self.tmin, self.tmax, nt).tolist()
826 @property
827 def deltat(self):
828 return self.deltatdays * 24 * 3600
831class PsCmpConfig(Object):
833 version = String.T(default='2008a')
834 # scatter, profile or array
835 observation = PsCmpObservation.T(default=PsCmpScatter.D())
837 rectangular_fault_size_factor = Float.T(
838 default=1.,
839 help='The size of the rectangular faults in the compound MT'
840 ' :py:class:`PsCmpMomentTensor` is dependend on the horizontal'
841 ' spacing of the GF database. This factor is applied to the'
842 ' relationship in i.e. fault length & width = f * dx.')
844 snapshots = PsCmpSnapshots.T(
845 optional=True)
847 rectangular_source_patches = List.T(PsCmpRectangularSource.T())
849 def items(self):
850 return dict(self.T.inamevals(self))
853class PsCmpConfigFull(PsCmpConfig):
855 pscmp_outdir = String.T(default='./')
856 psgrn_outdir = String.T(default='./psgrn_functions/')
858 los_vector = Tuple.T(3, Float.T(), optional=True)
860 sw_los_displacement = Int.T(default=0)
861 sw_coulomb_stress = Int.T(default=0)
862 coulomb_master_field = PsCmpCoulombStress.T(
863 optional=True,
864 default=PsCmpCoulombStressMasterFault.D())
866 displ_sw_output_types = Tuple.T(3, Int.T(), default=(0, 0, 0))
867 stress_sw_output_types = Tuple.T(6, Int.T(), default=(0, 0, 0, 0, 0, 0))
868 tilt_sw_output_types = Tuple.T(3, Int.T(), default=(0, 0, 0))
869 gravity_sw_output_types = Tuple.T(2, Int.T(), default=(0, 0))
871 indispl_filenames = Tuple.T(3, String.T(), default=psgrn_displ_names)
872 instress_filenames = Tuple.T(6, String.T(), default=psgrn_stress_names)
873 intilt_filenames = Tuple.T(3, String.T(), default=psgrn_tilt_names)
874 ingravity_filenames = Tuple.T(2, String.T(), default=psgrn_gravity_names)
876 outdispl_filenames = Tuple.T(3, String.T(), default=pscmp_displ_names)
877 outstress_filenames = Tuple.T(6, String.T(), default=pscmp_stress_names)
878 outtilt_filenames = Tuple.T(3, String.T(), default=pscmp_tilt_names)
879 outgravity_filenames = Tuple.T(2, String.T(), default=pscmp_gravity_names)
881 snapshot_basefilename = String.T(default='snapshot')
883 @classmethod
884 def example(cls):
885 conf = cls()
886 conf.psgrn_outdir = 'TEST_psgrn_functions/'
887 conf.pscmp_outdir = 'TEST_pscmp_output/'
888 conf.rectangular_source_patches = [PsCmpRectangularSource(
889 lat=10., lon=10., slip=2.,
890 width=5., length=10.,
891 strike=45, dip=30, rake=-90)]
892 conf.observation = PsCmpArray(
893 slat=9.5, elat=10.5, n_steps_lat=150,
894 slon=9.5, elon=10.5, n_steps_lon=150)
895 return conf
897 def get_output_filenames(self, rundir):
898 return [pjoin(rundir,
899 self.snapshot_basefilename + '_' + str(i + 1) + '.txt')
900 for i in range(len(self.snapshots.times))]
902 def string_for_config(self):
903 d = self.__dict__.copy()
905 patches_str, n_patches = distributed_fault_patches_to_config(
906 self.rectangular_source_patches)
908 d['patches_str'] = patches_str
909 d['n_patches'] = n_patches
911 str_npoints, str_observation = self.observation.string_for_config()
912 d['str_npoints'] = str_npoints
913 d['str_observation'] = str_observation
914 d['sw_observation_type'] = self.observation.sw
916 if self.snapshots.times:
917 str_times_dummy = []
918 for i, time in enumerate(self.snapshots.times):
919 str_times_dummy.append("%f '%s_%i.txt'\n" % (
920 time, self.snapshot_basefilename, i + 1))
922 str_times_dummy.append('#')
923 d['str_times_snapshots'] = ''.join(str_times_dummy)
924 d['n_snapshots'] = len(str_times_dummy) - 1
925 else:
926 d['str_times_snapshots'] = '# '
927 d['n_snapshots'] = 0
929 if self.sw_los_displacement == 1:
930 d['str_los_vector'] = str_float_vals(self.los_vector)
931 else:
932 d['str_los_vector'] = ''
934 if self.sw_coulomb_stress == 1:
935 d['str_coulomb_master_field'] = \
936 self.coulomb_master_field.string_for_config()
937 else:
938 d['str_coulomb_master_field'] = ''
940 d['str_psgrn_outdir'] = "'%s'" % self.psgrn_outdir
941 d['str_pscmp_outdir'] = "'%s'" % './'
943 d['str_indispl_filenames'] = str_str_vals(self.indispl_filenames)
944 d['str_instress_filenames'] = str_str_vals(self.instress_filenames)
945 d['str_intilt_filenames'] = str_str_vals(self.intilt_filenames)
946 d['str_ingravity_filenames'] = str_str_vals(self.ingravity_filenames)
948 d['str_outdispl_filenames'] = str_str_vals(self.outdispl_filenames)
949 d['str_outstress_filenames'] = str_str_vals(self.outstress_filenames)
950 d['str_outtilt_filenames'] = str_str_vals(self.outtilt_filenames)
951 d['str_outgravity_filenames'] = str_str_vals(self.outgravity_filenames)
953 d['str_displ_sw_output_types'] = str_int_vals(
954 self.displ_sw_output_types)
955 d['str_stress_sw_output_types'] = str_int_vals(
956 self.stress_sw_output_types)
957 d['str_tilt_sw_output_types'] = str_int_vals(
958 self.tilt_sw_output_types)
959 d['str_gravity_sw_output_types'] = str_int_vals(
960 self.gravity_sw_output_types)
962 template = '''# autogenerated PSCMP input by pscmp.py
963#===============================================================================
964# This is input file of FORTRAN77 program "pscmp08" for modeling post-seismic
965# deformation induced by earthquakes in multi-layered viscoelastic media using
966# the Green's function approach. The earthquke source is represented by an
967# arbitrary number of rectangular dislocation planes. For more details, please
968# read the accompanying READ.ME file.
969#
970# written by Rongjiang Wang
971# GeoForschungsZentrum Potsdam
972# e-mail: wang@gfz-potsdam.de
973# phone +49 331 2881209
974# fax +49 331 2881204
975#
976# Last modified: Potsdam, July, 2008
977#
978# References:
979#
980# (1) Wang, R., F. Lorenzo-Martin and F. Roth (2003), Computation of deformation
981# induced by earthquakes in a multi-layered elastic crust - FORTRAN programs
982# EDGRN/EDCMP, Computer and Geosciences, 29(2), 195-207.
983# (2) Wang, R., F. Lorenzo-Martin and F. Roth (2006), PSGRN/PSCMP - a new code for
984# calculating co- and post-seismic deformation, geoid and gravity changes
985# based on the viscoelastic-gravitational dislocation theory, Computers and
986# Geosciences, 32, 527-541. DOI:10.1016/j.cageo.2005.08.006.
987# (3) Wang, R. (2005), The dislocation theory: a consistent way for including the
988# gravity effect in (visco)elastic plane-earth models, Geophysical Journal
989# International, 161, 191-196.
990#
991#################################################################
992## ##
993## Green's functions should have been prepared with the ##
994## program "psgrn08" before the program "pscmp08" is started. ##
995## ##
996## For local Cartesian coordinate system, the Aki's convention ##
997## is used, that is, x is northward, y is eastward, and z is ##
998## downward. ##
999## ##
1000## If not specified otherwise, SI Unit System is used overall! ##
1001## ##
1002#################################################################
1003#===============================================================================
1004# OBSERVATION ARRAY
1005# =================
1006# 1. selection for irregular observation positions (= 0) or a 1D observation
1007# profile (= 1) or a rectangular 2D observation array (= 2): iposrec
1008#
1009# IF (iposrec = 0 for irregular observation positions) THEN
1010#
1011# 2. number of positions: nrec
1012#
1013# 3. coordinates of the observations: (lat(i),lon(i)), i=1,nrec
1014#
1015# ELSE IF (iposrec = 1 for regular 1D observation array) THEN
1016#
1017# 2. number of position samples of the profile: nrec
1018#
1019# 3. the start and end positions: (lat1,lon1), (lat2,lon2)
1020#
1021# ELSE IF (iposrec = 2 for rectanglular 2D observation array) THEN
1022#
1023# 2. number of x samples, start and end values: nxrec, xrec1, xrec2
1024#
1025# 3. number of y samples, start and end values: nyrec, yrec1, yrec2
1026#
1027# sequence of the positions in output data: lat(1),lon(1); ...; lat(nx),lon(1);
1028# lat(1),lon(2); ...; lat(nx),lon(2); ...; lat(1),lon(ny); ...; lat(nx),lon(ny).
1029#
1030# Note that the total number of observation positions (nrec or nxrec*nyrec)
1031# should be <= NRECMAX (see pecglob.h)!
1032#===============================================================================
1033 %(sw_observation_type)i
1034%(str_npoints)s
1035%(str_observation)s
1036#===============================================================================
1037# OUTPUTS
1038# =======
1039#
1040# 1. select output for los displacement (only for snapshots, see below), x, y,
1041# and z-cosines to the INSAR orbit: insar (1/0 = yes/no), xlos, ylos, zlos
1042#
1043# if this option is selected (insar = 1), the snapshots will include additional
1044# data:
1045# LOS_Dsp = los displacement to the given satellite orbit.
1046#
1047# 2. select output for Coulomb stress changes (only for snapshots, see below):
1048# icmb (1/0 = yes/no), friction, Skempton ratio, strike, dip, and rake angles
1049# [deg] describing the uniform regional master fault mechanism, the uniform
1050# regional principal stresses: sigma1, sigma2 and sigma3 [Pa] in arbitrary
1051# order (the orietation of the pre-stress field will be derived by assuming
1052# that the master fault is optimally oriented according to Coulomb failure
1053# criterion)
1054#
1055# if this option is selected (icmb = 1), the snapshots will include additional
1056# data:
1057# CMB_Fix, Sig_Fix = Coulomb and normal stress changes on master fault;
1058# CMB_Op1/2, Sig_Op1/2 = Coulomb and normal stress changes on the two optimally
1059# oriented faults;
1060# Str_Op1/2, Dip_Op1/2, Slp_Op1/2 = strike, dip and rake angles of the two
1061# optimally oriented faults.
1062#
1063# Note: the 1. optimally orieted fault is the one closest to the master fault.
1064#
1065# 3. output directory in char format: outdir
1066#
1067# 4. select outputs for displacement components (1/0 = yes/no): itout(i), i=1-3
1068#
1069# 5. the file names in char format for the x, y, and z components:
1070# toutfile(i), i=1-3
1071#
1072# 6. select outputs for stress components (1/0 = yes/no): itout(i), i=4-9
1073#
1074# 7. the file names in char format for the xx, yy, zz, xy, yz, and zx components:
1075# toutfile(i), i=4-9
1076#
1077# 8. select outputs for vertical NS and EW tilt components, block rotation, geoid
1078# and gravity changes (1/0 = yes/no): itout(i), i=10-14
1079#
1080# 9. the file names in char format for the NS tilt (positive if borehole top
1081# tilts to north), EW tilt (positive if borehole top tilts to east), block
1082# rotation (clockwise positive), geoid and gravity changes: toutfile(i), i=10-14
1083#
1084# Note that all above outputs are time series with the time window as same
1085# as used for the Green's functions
1086#
1087#10. number of scenario outputs ("snapshots": spatial distribution of all above
1088# observables at given time points; <= NSCENMAX (see pscglob.h): nsc
1089#
1090#11. the time [day], and file name (in char format) for the 1. snapshot;
1091#12. the time [day], and file name (in char format) for the 2. snapshot;
1092#13. ...
1093#
1094# Note that all file or directory names should not be longer than 80
1095# characters. Directories must be ended by / (unix) or \ (dos)!
1096#===============================================================================
1097 %(sw_los_displacement)i %(str_los_vector)s
1098 %(sw_coulomb_stress)i %(str_coulomb_master_field)s
1099 %(str_pscmp_outdir)s
1100 %(str_displ_sw_output_types)s
1101 %(str_outdispl_filenames)s
1102 %(str_stress_sw_output_types)s
1103 %(str_outstress_filenames)s
1104 %(str_tilt_sw_output_types)s %(str_gravity_sw_output_types)s
1105 %(str_outtilt_filenames)s %(str_outgravity_filenames)s
1106 %(n_snapshots)i
1107%(str_times_snapshots)s
1108#===============================================================================
1109#
1110# GREEN'S FUNCTION DATABASE
1111# =========================
1112# 1. directory where the Green's functions are stored: grndir
1113#
1114# 2. file names (without extensions!) for the 13 Green's functions:
1115# 3 displacement komponents (uz, ur, ut): green(i), i=1-3
1116# 6 stress components (szz, srr, stt, szr, srt, stz): green(i), i=4-9
1117# radial and tangential components measured by a borehole tiltmeter,
1118# rigid rotation around z-axis, geoid and gravity changes (tr, tt, rot, gd, gr):
1119# green(i), i=10-14
1120#
1121# Note that all file or directory names should not be longer than 80
1122# characters. Directories must be ended by / (unix) or \ (dos)! The
1123# extensions of the file names will be automatically considered. They
1124# are ".ep", ".ss", ".ds" and ".cl" denoting the explosion (inflation)
1125# strike-slip, the dip-slip and the compensated linear vector dipole
1126# sources, respectively.
1127#
1128#===============================================================================
1129 %(str_psgrn_outdir)s
1130 %(str_indispl_filenames)s
1131 %(str_instress_filenames)s
1132 %(str_intilt_filenames)s %(str_ingravity_filenames)s
1133#===============================================================================
1134# RECTANGULAR SUBFAULTS
1135# =====================
1136# 1. number of subfaults (<= NSMAX in pscglob.h): ns
1137#
1138# 2. parameters for the 1. rectangular subfault: geographic coordinates
1139# (O_lat, O_lon) [deg] and O_depth [km] of the local reference point on
1140# the present fault plane, length (along strike) [km] and width (along down
1141# dip) [km], strike [deg], dip [deg], number of equi-size fault patches along
1142# the strike (np_st) and along the dip (np_di) (total number of fault patches
1143# = np_st x np_di), and the start time of the rupture; the following data
1144# lines describe the slip distribution on the present sub-fault:
1145#
1146# pos_s[km] pos_d[km] slip_strike[m] slip_downdip[m] opening[m]
1147#
1148# where (pos_s,pos_d) defines the position of the center of each patch in
1149# the local coordinate system with the origin at the reference point:
1150# pos_s = distance along the length (positive in the strike direction)
1151# pos_d = distance along the width (positive in the down-dip direction)
1152#
1153#
1154# 3. ... for the 2. subfault ...
1155# ...
1156# N
1157# /
1158# /| strike
1159# +------------------------
1160# |\ p . \ W
1161# :-\ i . \ i
1162# | \ l . \ d
1163# :90 \ S . \ t
1164# |-dip\ . \ h
1165# : \. | rake \
1166# Z -------------------------
1167# L e n g t h
1168#
1169# Simulation of a Mogi source:
1170# (1) Calculate deformation caused by three small openning plates (each
1171# causes a third part of the volume of the point inflation) located
1172# at the same depth as the Mogi source but oriented orthogonal to
1173# each other.
1174# (2) Multiply the results by 3(1-nu)/(1+nu), where nu is the Poisson
1175# ratio at the source depth.
1176# The multiplication factor is the ratio of the seismic moment (energy) of
1177# the Mogi source to that of the plate openning with the same volume change.
1178#===============================================================================
1179# n_faults
1180#-------------------------------------------------------------------------------
1181 %(n_patches)i
1182#-------------------------------------------------------------------------------
1183# n O_lat O_lon O_depth length width strike dip np_st np_di start_time
1184# [-] [deg] [deg] [km] [km] [km] [deg] [deg] [-] [-] [day]
1185# pos_s pos_d slp_stk slp_ddip open
1186# [km] [km] [m] [m] [m]
1187#-------------------------------------------------------------------------------
1188%(patches_str)s
1189#================================end of input===================================
1190''' # noqa
1191 return (template % d).encode('ascii')
1194class PsGrnPsCmpConfig(Object):
1195 '''
1196 Combined config Object of PsGrn and PsCmp.
1198 Note: attributes in this class use non-standard units [days] to be
1199 consistent with PSCMP text file input. Please read the documentation
1200 carefully.
1201 '''
1202 tmin_days = Float.T(
1203 default=0.,
1204 help='Min. time in days')
1205 tmax_days = Float.T(
1206 default=1.,
1207 help='Max. time in days')
1208 gf_outdir = String.T(default='psgrn_functions')
1210 psgrn_config = PsGrnConfig.T(default=PsGrnConfig.D())
1211 pscmp_config = PsCmpConfig.T(default=PsCmpConfig.D())
1214class PsCmpError(gf.store.StoreError):
1215 pass
1218class Interrupted(gf.store.StoreError):
1219 def __str__(self):
1220 return 'Interrupted.'
1223class PsCmpRunner(object):
1224 '''
1225 Wrapper object to execute the program fomosto_pscmp with the specified
1226 configuration.
1228 :param tmp: string, path to the temporary directy where calculation
1229 results are stored
1230 :param keep_tmp: boolean, if True the result directory is kept
1231 '''
1232 def __init__(self, tmp=None, keep_tmp=False):
1233 if tmp is not None:
1234 tmp = os.path.abspath(tmp)
1235 self.tempdir = mkdtemp(prefix='pscmprun-', dir=tmp)
1236 self.keep_tmp = keep_tmp
1237 self.config = None
1239 def run(self, config):
1240 '''
1241 Run the program!
1243 :param config: :py:class:`PsCmpConfigFull`
1244 '''
1245 self.config = config
1247 input_fn = pjoin(self.tempdir, 'input')
1249 with open(input_fn, 'wb') as f:
1250 input_str = config.string_for_config()
1252 logger.debug('===== begin pscmp input =====\n'
1253 '%s===== end pscmp input =====' % input_str.decode())
1255 f.write(input_str)
1257 program = program_bins['pscmp.%s' % config.version]
1259 old_wd = os.getcwd()
1260 os.chdir(self.tempdir)
1262 interrupted = []
1264 def signal_handler(signum, frame):
1265 os.kill(proc.pid, signal.SIGTERM)
1266 interrupted.append(True)
1268 original = signal.signal(signal.SIGINT, signal_handler)
1269 try:
1270 try:
1271 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE,
1272 close_fds=True)
1274 except OSError as err:
1275 os.chdir(old_wd)
1276 logger.error('OS error: {0}'.format(err))
1277 raise PsCmpError(
1278 '''could not start pscmp executable: "%s"
1279Available fomosto backends and download links to the modelling codes are listed
1280on
1282 https://pyrocko.org/docs/current/apps/fomosto/backends.html
1284''' % program)
1286 (output_str, error_str) = proc.communicate(b'input\n')
1288 finally:
1289 signal.signal(signal.SIGINT, original)
1291 if interrupted:
1292 raise KeyboardInterrupt()
1294 logger.debug('===== begin pscmp output =====\n'
1295 '%s===== end pscmp output =====' % output_str.decode())
1297 errmsg = []
1298 if proc.returncode != 0:
1299 errmsg.append(
1300 'pscmp had a non-zero exit state: %i' % proc.returncode)
1302 if error_str:
1303 errmsg.append('pscmp emitted something via stderr')
1305 if output_str.lower().find(b'error') != -1:
1306 errmsg.append("the string 'error' appeared in pscmp output")
1308 if errmsg:
1309 self.keep_tmp = True
1311 os.chdir(old_wd)
1312 raise PsCmpError('''
1313===== begin pscmp input =====
1314{pscmp_input}===== end pscmp input =====
1315===== begin pscmp output =====
1316{pscmp_output}===== end pscmp output =====
1317===== begin pscmp error =====
1318{pscmp_error}===== end pscmp error =====
1319{error_messages}
1320pscmp has been invoked as "{call}"
1321in the directory {dir}'''.format(
1322 pscmp_input=input_str,
1323 pscmp_output=output_str,
1324 pscmp_error=error_str,
1325 error_messages='\n'.join(errmsg),
1326 call=program,
1327 dir=self.tempdir)
1328 .lstrip())
1330 self.pscmp_output = output_str
1331 self.pscmp_error = error_str
1333 os.chdir(old_wd)
1335 def get_results(self, component='displ'):
1336 '''
1337 Get the resulting components from the stored directory.
1338 Be careful: The z-component is downward positive!
1340 :param component: string, the component to retrieve from the
1341 result directory, may be:
1342 "displ": displacement, n x 3 array
1343 "stress": stresses n x 6 array
1344 "tilt': tilts n x 3 array,
1345 "gravity': gravity n x 2 array
1346 "all": all the above together
1347 '''
1348 assert self.config.snapshots is not None
1349 fns = self.config.get_output_filenames(self.tempdir)
1351 output = []
1352 for fn in fns:
1353 if not os.path.exists(fn):
1354 continue
1356 data = num.loadtxt(fn, skiprows=1, dtype=float)
1358 try:
1359 _, idxs = pscmp_component_mapping[component]
1360 except KeyError:
1361 raise Exception('component %s not supported! Either: %s' % (
1362 component, ', '.join(
1363 '"%s"' % k for k in pscmp_component_mapping.keys())))
1365 istart, iend = idxs
1367 output.append(data[:, istart:iend])
1369 return output
1371 def get_traces(self, component='displ'):
1372 '''
1373 Load snapshot data array and return specified components.
1374 Transform array component and receiver wise to list of
1375 :py:class:`pyrocko.trace.Trace`
1376 '''
1378 distances = self.config.observation.distances
1380 output_list = self.get_results(component=component)
1381 deltat = self.config.snapshots.deltat
1383 nrec, ncomp = output_list[0].shape
1385 outarray = num.dstack([num.zeros([nrec, ncomp])] + output_list)
1387 comps, _ = pscmp_component_mapping[component]
1389 outtraces = []
1390 for row in range(nrec):
1391 for col, comp in enumerate(comps):
1392 tr = trace.Trace(
1393 '', '%04i' % row, '', comp,
1394 tmin=0.0, deltat=deltat, ydata=outarray[row, col, :],
1395 meta=dict(distance=distances[row]))
1396 outtraces.append(tr)
1398 return outtraces
1400 def __del__(self):
1401 if self.tempdir:
1402 if not self.keep_tmp:
1403 shutil.rmtree(self.tempdir)
1404 self.tempdir = None
1405 else:
1406 logger.warning(
1407 'not removing temporary directory: %s' % self.tempdir)
1410class PsGrnCmpGFBuilder(gf.builder.Builder):
1411 nsteps = 2
1413 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
1414 force=False):
1416 self.store = gf.store.Store(store_dir, 'w')
1418 storeconf = self.store.config
1420 dummy_lat = 10.0
1421 dummy_lon = 10.0
1423 depths = storeconf.coords[0]
1424 lats = num.ones_like(depths) * dummy_lat
1425 lons = num.ones_like(depths) * dummy_lon
1426 points = num.vstack([lats, lons, depths]).T
1428 self.shear_moduli = storeconf.get_shear_moduli(
1429 lat=dummy_lat,
1430 lon=dummy_lon,
1431 points=points,
1432 interpolation='multilinear')
1434 self.lambda_moduli = storeconf.get_lambda_moduli(
1435 lat=dummy_lat,
1436 lon=dummy_lon,
1437 points=points,
1438 interpolation='multilinear')
1440 if step == 0:
1441 block_size = (1, storeconf.nsource_depths, storeconf.ndistances)
1442 else:
1443 if block_size is None:
1444 block_size = (1, 1, storeconf.ndistances)
1446 if len(storeconf.ns) == 2:
1447 block_size = block_size[1:]
1449 gf.builder.Builder.__init__(
1450 self, storeconf, step, block_size=block_size, force=force)
1452 baseconf = self.store.get_extra('psgrn_pscmp')
1454 cg = PsGrnConfigFull(**baseconf.psgrn_config.items())
1455 if cg.n_time_samples is None or cg.max_time is None:
1456 deltat_days = 1. / storeconf.sample_rate / day
1458 cg.n_time_samples = int(baseconf.tmax_days // deltat_days)
1459 cg.max_time = baseconf.tmax_days
1460 else:
1461 warnings.warn(
1462 'PsGrnConfig is defining n_times_samples and max_time,'
1463 ' this is replaced by PsGrnPsCmpConfig tmin and tmax.',
1464 FutureWarning)
1466 cc = PsCmpConfigFull(**baseconf.pscmp_config.items())
1467 if cc.snapshots is None:
1468 deltat_days = 1. / storeconf.sample_rate / day
1470 cc.snapshots = PsCmpSnapshots(
1471 tmin=baseconf.tmin_days,
1472 tmax=baseconf.tmax_days,
1473 deltatdays=deltat_days)
1474 else:
1475 warnings.warn(
1476 'PsCmpConfig is defining snapshots,'
1477 ' this is replaced by PsGrnPsCmpConfig tmin and tmax.',
1478 FutureWarning)
1480 cg.earthmodel_1d = storeconf.earthmodel_1d
1482 gf_outpath = os.path.join(store_dir, baseconf.gf_outdir)
1484 cg.psgrn_outdir = gf_outpath + '/'
1485 cc.psgrn_outdir = gf_outpath + '/'
1487 util.ensuredir(gf_outpath)
1489 self.cg = cg
1490 self.cc = cc
1492 def cleanup(self):
1493 self.store.close()
1495 def work_block(self, iblock):
1497 if len(self.store.config.ns) == 2:
1498 (sz, firstx), (sz, lastx), (ns, nx) = \
1499 self.get_block_extents(iblock)
1500 mu = self.shear_moduli[iblock]
1501 lame_lambda = self.lambda_moduli[iblock]
1503 rz = self.store.config.receiver_depth
1504 else:
1505 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
1506 self.get_block_extents(iblock)
1508 fc = self.store.config
1509 cc = copy.deepcopy(self.cc)
1510 cg = copy.deepcopy(self.cg)
1512 dx = fc.distance_delta
1514 logger.info(
1515 'Starting step %i / %i, block %i / %i' %
1516 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
1518 if self.step == 0:
1519 if cg.gf_depth_spacing is None:
1520 gf_depth_spacing = fc.source_depth_delta
1521 else:
1522 gf_depth_spacing = cg.gf_depth_spacing * km
1524 n_steps_depth = int((fc.source_depth_max - fc.source_depth_min) /
1525 gf_depth_spacing) + 1
1527 if cg.gf_distance_spacing is None:
1528 gf_distance_spacing = fc.distance_delta
1529 else:
1530 gf_distance_spacing = cg.gf_distance_spacing * km
1531 n_steps_distance = int((fc.distance_max - fc.distance_min) /
1532 gf_distance_spacing) + 1
1534 cg.depth_grid = PsGrnSpatialSampling(
1535 n_steps=n_steps_depth,
1536 start_distance=fc.source_depth_min / km,
1537 end_distance=fc.source_depth_max / km)
1539 cg.distance_grid = PsGrnSpatialSampling(
1540 n_steps=n_steps_distance,
1541 start_distance=fc.distance_min / km,
1542 end_distance=fc.distance_max / km)
1544 runner = PsGrnRunner(outdir=self.cg.psgrn_outdir)
1545 runner.run(cg, force=self.force)
1547 else:
1548 distances = num.linspace(
1549 firstx, firstx + (nx - 1) * dx, nx).tolist()
1551 # fomosto sample rate in s, pscmp takes days
1552 deltatdays = 1. / (fc.sample_rate * 24. * 3600.)
1553 cc.snapshots = PsCmpSnapshots(
1554 tmin=0., tmax=cg.max_time, deltatdays=deltatdays)
1555 cc.observation = PsCmpProfile(
1556 slat=0. - 0.001 * cake.m2d,
1557 slon=0.0,
1558 elat=0. + distances[-1] * cake.m2d,
1559 elon=0.0,
1560 n_steps=len(distances),
1561 distances=distances)
1563 runner = PsCmpRunner(keep_tmp=False)
1565 mtsize = float(dx * cc.rectangular_fault_size_factor)
1567 Ai = 1. / num.power(mtsize, 2)
1568 nullf, sf = get_iso_scaling_factors(mu=mu, lame_lambda=lame_lambda)
1570 mui = 1. / mu
1572 gfmapping = [
1573 (('nn',),
1574 {'un': (0, Ai * sf), 'ud': (5, Ai * sf)}),
1575 (('ne',),
1576 {'ue': (3, 1 * Ai * mui)}),
1577 (('nd',),
1578 {'un': (1, 1 * Ai * mui), 'ud': (6, 1 * Ai * mui)}),
1579 (('ed',),
1580 {'ue': (4, 1 * Ai * mui)}),
1581 (('dd',),
1582 {'un': (2, Ai * sf), 'ud': (7, Ai * sf)}),
1583 (('ee',),
1584 {'un': (8, Ai * sf), 'ud': (9, Ai * sf)}),
1585 ]
1587 for mt, gfmap in gfmapping:
1588 cc.rectangular_source_patches = []
1589 for idx in mt:
1590 pmt = PsCmpMomentTensor(
1591 lat=0. + 0.001 * dx * cake.m2d,
1592 lon=0.0,
1593 depth=float(sz),
1594 width=mtsize,
1595 length=mtsize,
1596 idx=idx)
1598 cc.rectangular_source_patches.extend(
1599 pmt.to_rfs(nullf))
1601 runner.run(cc)
1603 rawtraces = runner.get_traces()
1605 interrupted = []
1607 def signal_handler(signum, frame):
1608 interrupted.append(True)
1610 original = signal.signal(signal.SIGINT, signal_handler)
1611 self.store.lock()
1612 duplicate_inserts = 0
1614 try:
1615 for itr, tr in enumerate(rawtraces):
1616 if tr.channel in gfmap:
1617 x = tr.meta['distance']
1618 ig, factor = gfmap[tr.channel]
1620 if len(self.store.config.ns) == 2:
1621 args = (sz, x, ig)
1622 else:
1623 args = (rz, sz, x, ig)
1625 gf_tr = gf.store.GFTrace.from_trace(tr)
1627 gf_tr.data *= factor
1629 try:
1630 self.store.put(args, gf_tr)
1631 except gf.store.DuplicateInsert:
1632 duplicate_inserts += 1
1634 finally:
1635 if duplicate_inserts:
1636 logger.warning(
1637 '%i insertions skipped (duplicates)'
1638 % duplicate_inserts)
1640 self.store.unlock()
1641 signal.signal(signal.SIGINT, original)
1643 if interrupted:
1644 raise KeyboardInterrupt()
1646 logger.info(
1647 'Done with step %i / %i, block %i / %i' % (
1648 self.step + 1, self.nsteps, iblock + 1, self.nblocks))
1651def init(store_dir, variant):
1652 if variant is None:
1653 variant = '2008a'
1655 if ('pscmp.' + variant) not in program_bins:
1656 raise gf.store.StoreError('unsupported pscmp variant: %s' % variant)
1658 if ('psgrn.' + variant) not in program_bins:
1659 raise gf.store.StoreError('unsupported psgrn variant: %s' % variant)
1661 c = PsGrnPsCmpConfig()
1663 store_id = os.path.basename(os.path.realpath(store_dir))
1665 # Initialising a viscous mantle
1666 cake_mod = cake.load_model(fn=None, crust2_profile=(54., 23.))
1667 mantle = cake_mod.material(z=45*km)
1668 mantle.burger_eta1 = 5e17
1669 mantle.burger_eta2 = 1e19
1670 mantle.burger_alpha = 1.
1672 config = gf.meta.ConfigTypeA(
1673 id=store_id,
1674 ncomponents=10,
1675 sample_rate=1. / (3600. * 24.),
1676 receiver_depth=0. * km,
1677 source_depth_min=0. * km,
1678 source_depth_max=15. * km,
1679 source_depth_delta=.5 * km,
1680 distance_min=0. * km,
1681 distance_max=50. * km,
1682 distance_delta=1. * km,
1683 earthmodel_1d=cake_mod,
1684 modelling_code_id='psgrn_pscmp.%s' % variant,
1685 tabulated_phases=[]) # dummy list
1687 c.validate()
1688 config.validate()
1689 return gf.store.Store.create_editables(
1690 store_dir,
1691 config=config,
1692 extra={'psgrn_pscmp': c})
1695def build(store_dir,
1696 force=False,
1697 nworkers=None,
1698 continue_=False,
1699 step=None,
1700 iblock=None):
1702 return PsGrnCmpGFBuilder.build(
1703 store_dir, force=force, nworkers=nworkers, continue_=continue_,
1704 step=step, iblock=iblock)