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-04 09:52 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import logging
7import warnings
8import os
9from os.path import join as pjoin
10import signal
11import shutil
12import copy
14import math
15import numpy as num
17from tempfile import mkdtemp
18from subprocess import Popen, PIPE
20from pyrocko.guts import Float, Int, Tuple, List, Object, String
21from pyrocko.model import Location
22from pyrocko import gf, util, trace, cake
25# how to call the programs
26program_bins = {
27 'pscmp.2008a': 'fomosto_pscmp2008a',
28 'psgrn.2008a': 'fomosto_psgrn2008a'
29}
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()
37km = 1000.
38day = 3600. * 24
40guts_prefix = 'pf'
41logger = logging.getLogger('pyrocko.fomosto.psgrn_pscmp')
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()
50 except OSError:
51 return False
53 return True
56def nextpow2(i):
57 return 2 ** int(math.ceil(math.log(i) / math.log(2.)))
60def str_float_vals(vals):
61 return ' '.join('%e' % val for val in vals)
64def str_int_vals(vals):
65 return ' '.join('%i' % val for val in vals)
68def str_str_vals(vals):
69 return ' '.join("'%s'" % val for val in vals)
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)))
80 return '\n'.join(srows), len(srows)
83class PsGrnSpatialSampling(Object):
85 '''
86 Definition of spatial sampling for PSGRN.
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 '''
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
96 def string_for_config(self):
97 return '%i %15e %15e' % (self.n_steps, self.start_distance,
98 self.end_distance)
101class PsGrnConfig(Object):
103 '''
104 Configuration for PSGRN.
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 '''
111 version = String.T(default='2008a')
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.')
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]')
137 def items(self):
138 return dict(self.T.inamevals(self))
141class PsGrnConfigFull(PsGrnConfig):
143 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True)
144 psgrn_outdir = String.T(default='psgrn_green/')
146 distance_grid = PsGrnSpatialSampling.T(PsGrnSpatialSampling.D())
147 depth_grid = PsGrnSpatialSampling.T(PsGrnSpatialSampling.D())
149 sw_source_regime = Int.T(default=1) # 1-continental, 0-ocean
150 sw_gravity = Int.T(default=0)
152 accuracy_wavenumber_integration = Float.T(default=0.025)
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)
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
166 def string_for_config(self):
168 assert self.earthmodel_1d is not None
170 d = self.__dict__.copy()
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
177 d['str_psgrn_outdir'] = "'%s'" % './'
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)
184 d['str_distance_grid'] = self.distance_grid.string_for_config()
185 d['str_depth_grid'] = self.depth_grid.string_for_config()
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')
329class PsGrnError(gf.store.StoreError):
330 pass
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)
341class PsGrnRunner(object):
342 '''
343 Wrapper object to execute the program fomosto_psgrn.
344 '''
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
353 def run(self, config, force=False):
354 '''
355 Run the program with the specified configuration.
357 :param config: :py:class:`PsGrnConfigFull`
358 :param force: boolean, set true to overwrite existing output
359 '''
360 self.config = config
362 input_fn = pjoin(self.outdir, 'input')
364 remove_if_exists(input_fn, force=force)
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]
373 old_wd = os.getcwd()
375 os.chdir(self.outdir)
377 interrupted = []
379 def signal_handler(signum, frame):
380 os.kill(proc.pid, signal.SIGTERM)
381 interrupted.append(True)
383 original = signal.signal(signal.SIGINT, signal_handler)
384 try:
385 try:
386 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
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
395 https://pyrocko.org/docs/current/apps/fomosto/backends.html
397''' % program)
399 (output_str, error_str) = proc.communicate(b'input\n')
401 finally:
402 signal.signal(signal.SIGINT, original)
404 if interrupted:
405 raise KeyboardInterrupt()
407 logger.debug('===== begin psgrn output =====\n'
408 '%s===== end psgrn output =====' % output_str.decode())
410 errmess = []
411 if proc.returncode != 0:
412 errmess.append(
413 'psgrn had a non-zero exit state: %i' % proc.returncode)
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')
421 if output_str.lower().find(b'error') != -1:
422 errmess.append("the string 'error' appeared in psgrn output")
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))
439 self.psgrn_output = output_str
440 self.psgrn_error = error_str
442 os.chdir(old_wd)
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
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 }
461def dsin(value):
462 return num.sin(value * num.pi / 180.)
465def dcos(value):
466 return num.cos(value * num.pi / 180.)
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()))
477 return '\n'.join(srows), len(srows)
480class PsCmpObservation(Object):
481 pass
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])
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))
496 self.sw = 0
497 return ' %i' % (len(srows)), '\n'.join(srows)
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.')
522 def string_for_config(self):
523 self.sw = 1
525 return ' %i' % (self.n_steps), \
526 ' ( %15f, %15f ), ( %15f, %15f )' % (
527 self.slat, self.slon, self.elat, self.elon)
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')
549 def string_for_config(self):
550 self.sw = 2
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)
558class PsCmpRectangularSource(Location, gf.seismosizer.Cloneable):
559 '''
560 Rectangular Source for the input geometry of the active fault.
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)
581 slip = Float.T(optional=True, default=1.0)
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)
587 def update(self, **kwargs):
588 '''
589 Change some of the source models parameters.
591 Example::
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
605 '''
606 for (k, v) in kwargs.items():
607 self[k] = v
609 @property
610 def dip_slip(self):
611 return float(self.slip * dsin(self.rake) * (-1))
613 @property
614 def strike_slip(self):
615 return float(self.slip * dcos(self.rake))
617 def string_for_config(self):
619 if self.pos_s or self.pos_d is None:
620 self.pos_s = 0.
621 self.pos_d = 0.
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
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
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 }
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 }
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)
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))
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
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!
677 Warning: Mixed standard [m] / non-standard [days] units are used in this
678 class. Please read the documentation carefully.
679 '''
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"')
694 def to_rfs(self, nullification):
696 cmpd = []
697 for comp, mt in MTIso.items():
698 params = copy.deepcopy(mt)
700 if comp != self.idx:
701 params = copy.deepcopy(mt)
702 params['opening'] *= nullification
704 kwargs = copy.deepcopy(self.__dict__)
705 kwargs.update(params)
706 kwargs.pop('idx')
707 kwargs.pop('_latlon')
708 cmpd.append(PsCmpRectangularSource(**kwargs))
710 return cmpd
713class PsCmpShearSF(Location, gf.seismosizer.Cloneable):
714 '''
715 Shear fault source model component.
717 Warning: Mixed standard [m] / non-standard [days] units are used in this
718 class. Please read the documentation carefully.
719 '''
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"')
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)]
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!
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".')
764 def to_rfs(self, nullification=-0.25):
765 kwargs = copy.deepcopy(self.__dict__)
766 kwargs.pop('_latlon')
768 if self.idx in MTIso:
769 tsf = PsCmpTensileSF(**kwargs)
770 return tsf.to_rfs(nullification)
772 elif self.idx in MTDev:
773 ssf = PsCmpShearSF(**kwargs)
774 return ssf.to_rfs()
776 else:
777 raise Exception('MT component not supported!')
780class PsCmpCoulombStress(Object):
781 pass
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]')
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__
800class PsCmpSnapshots(Object):
801 '''
802 Snapshot time series definition.
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 '''
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].')
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()
825 @property
826 def deltat(self):
827 return self.deltatdays * 24 * 3600
830class PsCmpConfig(Object):
832 version = String.T(default='2008a')
833 # scatter, profile or array
834 observation = PsCmpObservation.T(default=PsCmpScatter.D())
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.')
843 snapshots = PsCmpSnapshots.T(
844 optional=True)
846 rectangular_source_patches = List.T(PsCmpRectangularSource.T())
848 def items(self):
849 return dict(self.T.inamevals(self))
852class PsCmpConfigFull(PsCmpConfig):
854 pscmp_outdir = String.T(default='./')
855 psgrn_outdir = String.T(default='./psgrn_functions/')
857 los_vector = Tuple.T(3, Float.T(), optional=True)
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())
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))
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)
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)
880 snapshot_basefilename = String.T(default='snapshot')
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
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))]
901 def string_for_config(self):
902 d = self.__dict__.copy()
904 patches_str, n_patches = distributed_fault_patches_to_config(
905 self.rectangular_source_patches)
907 d['patches_str'] = patches_str
908 d['n_patches'] = n_patches
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
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))
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
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'] = ''
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'] = ''
939 d['str_psgrn_outdir'] = "'%s'" % self.psgrn_outdir
940 d['str_pscmp_outdir'] = "'%s'" % './'
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)
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)
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)
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')
1193class PsGrnPsCmpConfig(Object):
1194 '''
1195 Combined config Object of PsGrn and PsCmp.
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')
1209 psgrn_config = PsGrnConfig.T(default=PsGrnConfig.D())
1210 pscmp_config = PsCmpConfig.T(default=PsCmpConfig.D())
1213class PsCmpError(gf.store.StoreError):
1214 pass
1217class Interrupted(gf.store.StoreError):
1218 def __str__(self):
1219 return 'Interrupted.'
1222class PsCmpRunner(object):
1223 '''
1224 Wrapper object to execute the program fomosto_pscmp with the specified
1225 configuration.
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
1238 def run(self, config):
1239 '''
1240 Run the program!
1242 :param config: :py:class:`PsCmpConfigFull`
1243 '''
1244 self.config = config
1246 input_fn = pjoin(self.tempdir, 'input')
1248 with open(input_fn, 'wb') as f:
1249 input_str = config.string_for_config()
1251 logger.debug('===== begin pscmp input =====\n'
1252 '%s===== end pscmp input =====' % input_str.decode())
1254 f.write(input_str)
1256 program = program_bins['pscmp.%s' % config.version]
1258 old_wd = os.getcwd()
1259 os.chdir(self.tempdir)
1261 interrupted = []
1263 def signal_handler(signum, frame):
1264 os.kill(proc.pid, signal.SIGTERM)
1265 interrupted.append(True)
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)
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
1281 https://pyrocko.org/docs/current/apps/fomosto/backends.html
1283''' % program)
1285 (output_str, error_str) = proc.communicate(b'input\n')
1287 finally:
1288 signal.signal(signal.SIGINT, original)
1290 if interrupted:
1291 raise KeyboardInterrupt()
1293 logger.debug('===== begin pscmp output =====\n'
1294 '%s===== end pscmp output =====' % output_str.decode())
1296 errmsg = []
1297 if proc.returncode != 0:
1298 errmsg.append(
1299 'pscmp had a non-zero exit state: %i' % proc.returncode)
1301 if error_str:
1302 errmsg.append('pscmp emitted something via stderr')
1304 if output_str.lower().find(b'error') != -1:
1305 errmsg.append("the string 'error' appeared in pscmp output")
1307 if errmsg:
1308 self.keep_tmp = True
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())
1329 self.pscmp_output = output_str
1330 self.pscmp_error = error_str
1332 os.chdir(old_wd)
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!
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)
1350 output = []
1351 for fn in fns:
1352 if not os.path.exists(fn):
1353 continue
1355 data = num.loadtxt(fn, skiprows=1, dtype=float)
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())))
1364 istart, iend = idxs
1366 output.append(data[:, istart:iend])
1368 return output
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 '''
1377 distances = self.config.observation.distances
1379 output_list = self.get_results(component=component)
1380 deltat = self.config.snapshots.deltat
1382 nrec, ncomp = output_list[0].shape
1384 outarray = num.dstack([num.zeros([nrec, ncomp])] + output_list)
1386 comps, _ = pscmp_component_mapping[component]
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)
1397 return outtraces
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)
1409class PsGrnCmpGFBuilder(gf.builder.Builder):
1410 nsteps = 2
1412 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
1413 force=False):
1415 self.store = gf.store.Store(store_dir, 'w')
1417 storeconf = self.store.config
1419 dummy_lat = 10.0
1420 dummy_lon = 10.0
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
1427 self.shear_moduli = storeconf.get_shear_moduli(
1428 lat=dummy_lat,
1429 lon=dummy_lon,
1430 points=points,
1431 interpolation='multilinear')
1433 self.lambda_moduli = storeconf.get_lambda_moduli(
1434 lat=dummy_lat,
1435 lon=dummy_lon,
1436 points=points,
1437 interpolation='multilinear')
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)
1445 if len(storeconf.ns) == 2:
1446 block_size = block_size[1:]
1448 gf.builder.Builder.__init__(
1449 self, storeconf, step, block_size=block_size, force=force)
1451 baseconf = self.store.get_extra('psgrn_pscmp')
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
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)
1465 cc = PsCmpConfigFull(**baseconf.pscmp_config.items())
1466 if cc.snapshots is None:
1467 deltat_days = 1. / storeconf.sample_rate / day
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)
1479 cg.earthmodel_1d = storeconf.earthmodel_1d
1481 gf_outpath = os.path.join(store_dir, baseconf.gf_outdir)
1483 cg.psgrn_outdir = gf_outpath + '/'
1484 cc.psgrn_outdir = gf_outpath + '/'
1486 util.ensuredir(gf_outpath)
1488 self.cg = cg
1489 self.cc = cc
1491 def cleanup(self):
1492 self.store.close()
1494 def work_block(self, iblock):
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]
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)
1507 fc = self.store.config
1508 cc = copy.deepcopy(self.cc)
1509 cg = copy.deepcopy(self.cg)
1511 dx = fc.distance_delta
1513 logger.info(
1514 'Starting step %i / %i, block %i / %i' %
1515 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
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
1523 n_steps_depth = int((fc.source_depth_max - fc.source_depth_min) /
1524 gf_depth_spacing) + 1
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
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)
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)
1543 runner = PsGrnRunner(outdir=self.cg.psgrn_outdir)
1544 runner.run(cg, force=self.force)
1546 else:
1547 distances = num.linspace(
1548 firstx, firstx + (nx - 1) * dx, nx).tolist()
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)
1562 runner = PsCmpRunner(keep_tmp=False)
1564 mtsize = float(dx * cc.rectangular_fault_size_factor)
1566 Ai = 1. / num.power(mtsize, 2)
1567 nullf, sf = get_iso_scaling_factors(mu=mu, lame_lambda=lame_lambda)
1569 mui = 1. / mu
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 ]
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)
1597 cc.rectangular_source_patches.extend(
1598 pmt.to_rfs(nullf))
1600 runner.run(cc)
1602 rawtraces = runner.get_traces()
1604 interrupted = []
1606 def signal_handler(signum, frame):
1607 interrupted.append(True)
1609 original = signal.signal(signal.SIGINT, signal_handler)
1610 self.store.lock()
1611 duplicate_inserts = 0
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]
1619 if len(self.store.config.ns) == 2:
1620 args = (sz, x, ig)
1621 else:
1622 args = (rz, sz, x, ig)
1624 gf_tr = gf.store.GFTrace.from_trace(tr)
1626 gf_tr.data *= factor
1628 try:
1629 self.store.put(args, gf_tr)
1630 except gf.store.DuplicateInsert:
1631 duplicate_inserts += 1
1633 finally:
1634 if duplicate_inserts:
1635 logger.warning(
1636 '%i insertions skipped (duplicates)'
1637 % duplicate_inserts)
1639 self.store.unlock()
1640 signal.signal(signal.SIGINT, original)
1642 if interrupted:
1643 raise KeyboardInterrupt()
1645 logger.info(
1646 'Done with step %i / %i, block %i / %i' % (
1647 self.step + 1, self.nsteps, iblock + 1, self.nblocks))
1650def init(store_dir, variant):
1651 if variant is None:
1652 variant = '2008a'
1654 if ('pscmp.' + variant) not in program_bins:
1655 raise gf.store.StoreError('unsupported pscmp variant: %s' % variant)
1657 if ('psgrn.' + variant) not in program_bins:
1658 raise gf.store.StoreError('unsupported psgrn variant: %s' % variant)
1660 c = PsGrnPsCmpConfig()
1662 store_id = os.path.basename(os.path.realpath(store_dir))
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.
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
1686 c.validate()
1687 config.validate()
1688 return gf.store.Store.create_editables(
1689 store_dir,
1690 config=config,
1691 extra={'psgrn_pscmp': c})
1694def build(store_dir,
1695 force=False,
1696 nworkers=None,
1697 continue_=False,
1698 step=None,
1699 iblock=None):
1701 return PsGrnCmpGFBuilder.build(
1702 store_dir, force=force, nworkers=nworkers, continue_=continue_,
1703 step=step, iblock=iblock)