Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/fomosto/psgrn_pscmp.py: 61%
594 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-02-05 09:37 +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 The default shift of the origin (:gattr:`pos_s`, :gattr:`pos_d`) with
563 respect to the reference coordinates
564 (:gattr:`~pyrocko.model.location.Location.lat`,
565 :gattr:`~pyrocko.model.location.Location.lon`) is zero, which implies that
566 the reference is the center of the fault plane! The calculation point is
567 always the center of the fault-plane. Setting :gattr:`pos_s` or
568 :gattr:`pos_d` moves the fault point with respect to the origin along
569 strike and dip direction, respectively.
570 '''
572 length = Float.T(
573 default=6.0 * km,
574 help='Fault length [m].')
575 width = Float.T(
576 default=5.0 * km,
577 help='Fault width [m].')
578 strike = Float.T(
579 default=0.0,
580 help='Fault strike, clockwise from north [deg].')
581 dip = Float.T(
582 default=90.0,
583 help='Fault dip, downward from horizontal [deg].')
584 rake = Float.T(
585 default=0.0,
586 help='Rake angle, counter-clock wise from right-horizontal [deg].')
587 torigin = Float.T(
588 default=0.0,
589 help='Origin time [s].')
590 slip = Float.T(
591 optional=True,
592 default=1.0,
593 help='Slip amount [m]')
594 pos_s = Float.T(
595 optional=True,
596 default=None,
597 help='Origin offset along strike [m].')
598 pos_d = Float.T(
599 optional=True,
600 default=None,
601 help='Origin offset along dip [m].')
602 opening = Float.T(
603 default=0.0,
604 help='Fault opening [m].')
606 def update(self, **kwargs):
607 '''
608 Change some of the source models parameters.
610 Example::
612 >>> from pyrocko import gf
613 >>> s = gf.DCSource()
614 >>> s.update(strike=66., dip=33.)
615 >>> print(s)
616 --- !pf.DCSource
617 depth: 0.0
618 time: 1970-01-01 00:00:00
619 magnitude: 6.0
620 strike: 66.0
621 dip: 33.0
622 rake: 0.0
624 '''
625 for (k, v) in kwargs.items():
626 self[k] = v
628 @property
629 def dip_slip(self):
630 return float(self.slip * dsin(self.rake) * (-1))
632 @property
633 def strike_slip(self):
634 return float(self.slip * dcos(self.rake))
636 def string_for_config(self):
638 if self.pos_s or self.pos_d is None:
639 self.pos_s = 0.
640 self.pos_d = 0.
642 tempd = copy.deepcopy(self.__dict__)
643 tempd['dip_slip'] = self.dip_slip
644 tempd['strike_slip'] = self.strike_slip
645 tempd['effective_lat'] = self.effective_lat
646 tempd['effective_lon'] = self.effective_lon
647 tempd['depth'] /= km
648 tempd['length'] /= km
649 tempd['width'] /= km
651 return '%(effective_lat)15f %(effective_lon)15f %(depth)15f' \
652 '%(length)15f %(width)15f %(strike)15f' \
653 '%(dip)15f 1 1 %(torigin)15f \n %(pos_s)15f %(pos_d)15f ' \
654 '%(strike_slip)15f %(dip_slip)15f %(opening)15f' % tempd
657MTIso = {
658 'nn': dict(strike=90., dip=90., rake=0., slip=0., opening=1.),
659 'ee': dict(strike=0., dip=90., rake=0., slip=0., opening=1.),
660 'dd': dict(strike=0., dip=0., rake=-90., slip=0., opening=1.),
661 }
663MTDev = {
664 'ne': dict(strike=90., dip=90., rake=180., slip=1., opening=0.),
665 'nd': dict(strike=180., dip=0., rake=0., slip=1., opening=0.),
666 'ed': dict(strike=270., dip=0., rake=0., slip=1., opening=0.),
667 }
670def get_nullification_factor(mu, lame_lambda):
671 '''
672 Factor that needs to be multiplied to 2 of the tensile sources to
673 eliminate two of the isotropic components.
674 '''
675 return -lame_lambda / (2. * mu + 2. * lame_lambda)
678def get_trace_normalization_factor(mu, lame_lambda, nullification):
679 '''
680 Factor to be multiplied to elementary GF trace to have unit displacement.
681 '''
682 return 1. / ((2. * mu) + lame_lambda + (2. * lame_lambda * nullification))
685def get_iso_scaling_factors(mu, lame_lambda):
686 nullf = get_nullification_factor(mu, lame_lambda)
687 scale = get_trace_normalization_factor(mu, lame_lambda, nullf)
688 return nullf, scale
691class PsCmpTensileSF(Location, gf.seismosizer.Cloneable):
692 '''
693 Compound dislocation of 3 perpendicular, rectangular sources to approximate
694 an opening single force couple. NED coordinate system!
696 Warning: Mixed standard [m] / non-standard [days] units are used in this
697 class. Please read the documentation carefully.
698 '''
700 length = Float.T(
701 default=1.0 * km,
702 help='Length of source rectangle [m].')
703 width = Float.T(
704 default=1.0 * km,
705 help='Width of source rectangle [m].')
706 torigin = Float.T(
707 default=0.0,
708 help='Origin time [days]')
709 idx = String.T(
710 default='nn',
711 help='Axis index for opening direction; "nn", "ee", or "dd"')
713 def to_rfs(self, nullification):
715 cmpd = []
716 for comp, mt in MTIso.items():
717 params = copy.deepcopy(mt)
719 if comp != self.idx:
720 params = copy.deepcopy(mt)
721 params['opening'] *= nullification
723 kwargs = copy.deepcopy(self.__dict__)
724 kwargs.update(params)
725 kwargs.pop('idx')
726 kwargs.pop('_latlon')
727 cmpd.append(PsCmpRectangularSource(**kwargs))
729 return cmpd
732class PsCmpShearSF(Location, gf.seismosizer.Cloneable):
733 '''
734 Shear fault source model component.
736 Warning: Mixed standard [m] / non-standard [days] units are used in this
737 class. Please read the documentation carefully.
738 '''
740 length = Float.T(
741 default=1.0 * km,
742 help='Length of source rectangle [m].')
743 width = Float.T(
744 default=1.0 * km,
745 help='Width of source rectangle [m].')
746 torigin = Float.T(
747 default=0.0,
748 help='Origin time [days]')
749 idx = String.T(
750 default='ne',
751 help='Axis index for opening direction; "ne", "nd", or "ed"')
753 def to_rfs(self):
754 kwargs = copy.deepcopy(self.__dict__)
755 kwargs.pop('idx')
756 kwargs.pop('_latlon')
757 kwargs.update(MTDev[self.idx])
758 return [PsCmpRectangularSource(**kwargs)]
761class PsCmpMomentTensor(Location, gf.seismosizer.Cloneable):
762 '''
763 Mapping of Moment Tensor components to rectangular faults.
764 Only one component at a time valid! NED coordinate system!
766 Warning: Mixed standard [m] / non-standard [days] units are used in this
767 class. Please read the documentation carefully.
768 '''
769 length = Float.T(
770 default=1.0 * km,
771 help='Length of source rectangle [m].')
772 width = Float.T(
773 default=1.0 * km,
774 help='Width of source rectangle [m].')
775 torigin = Float.T(
776 default=0.0,
777 help='Origin time [days]')
778 idx = String.T(
779 default='nn',
780 help='Axis index for MT component; '
781 '"nn", "ee", "dd", "ne", "nd", or "ed".')
783 def to_rfs(self, nullification=-0.25):
784 kwargs = copy.deepcopy(self.__dict__)
785 kwargs.pop('_latlon')
787 if self.idx in MTIso:
788 tsf = PsCmpTensileSF(**kwargs)
789 return tsf.to_rfs(nullification)
791 elif self.idx in MTDev:
792 ssf = PsCmpShearSF(**kwargs)
793 return ssf.to_rfs()
795 else:
796 raise Exception('MT component not supported!')
799class PsCmpCoulombStress(Object):
800 pass
803class PsCmpCoulombStressMasterFault(PsCmpCoulombStress):
804 friction = Float.T(default=0.7)
805 skempton_ratio = Float.T(default=0.0)
806 master_fault_strike = Float.T(default=300.)
807 master_fault_dip = Float.T(default=15.)
808 master_fault_rake = Float.T(default=90.)
809 sigma1 = Float.T(default=1.e6, help='[Pa]')
810 sigma2 = Float.T(default=-1.e6, help='[Pa]')
811 sigma3 = Float.T(default=0.0, help='[Pa]')
813 def string_for_config(self):
814 return '%(friction)15e %(skempton_ratio)15e %(master_fault_strike)15f'\
815 '%(master_fault_dip)15f %(master_fault_rake)15f'\
816 '%(sigma1)15e %(sigma2)15e %(sigma3)15e' % self.__dict__
819class PsCmpSnapshots(Object):
820 '''
821 Snapshot time series definition.
823 Note: attributes in this class use non-standard units [days] to be
824 consistent with PSCMP text file input. Please read the documentation
825 carefully.
826 '''
828 tmin = Float.T(
829 default=0.0,
830 help='Time [days] after source time to start temporal sample'
831 ' snapshots.')
832 tmax = Float.T(
833 default=1.0,
834 help='Time [days] after source time to end temporal sample f.')
835 deltatdays = Float.T(
836 default=1.0,
837 help='Sample period [days].')
839 @property
840 def times(self):
841 nt = int(num.ceil((self.tmax - self.tmin) / self.deltatdays))
842 return num.linspace(self.tmin, self.tmax, nt).tolist()
844 @property
845 def deltat(self):
846 return self.deltatdays * 24 * 3600
849class PsCmpConfig(Object):
851 version = String.T(default='2008a')
852 # scatter, profile or array
853 observation = PsCmpObservation.T(default=PsCmpScatter.D())
855 rectangular_fault_size_factor = Float.T(
856 default=1.,
857 help='The size of the rectangular faults in the compound MT'
858 ' :py:class:`PsCmpMomentTensor` is dependend on the horizontal'
859 ' spacing of the GF database. This factor is applied to the'
860 ' relationship in i.e. fault length & width = f * dx.')
862 snapshots = PsCmpSnapshots.T(
863 optional=True)
865 rectangular_source_patches = List.T(PsCmpRectangularSource.T())
867 def items(self):
868 return dict(self.T.inamevals(self))
871class PsCmpConfigFull(PsCmpConfig):
873 pscmp_outdir = String.T(default='./')
874 psgrn_outdir = String.T(default='./psgrn_functions/')
876 los_vector = Tuple.T(3, Float.T(), optional=True)
878 sw_los_displacement = Int.T(default=0)
879 sw_coulomb_stress = Int.T(default=0)
880 coulomb_master_field = PsCmpCoulombStress.T(
881 optional=True,
882 default=PsCmpCoulombStressMasterFault.D())
884 displ_sw_output_types = Tuple.T(3, Int.T(), default=(0, 0, 0))
885 stress_sw_output_types = Tuple.T(6, Int.T(), default=(0, 0, 0, 0, 0, 0))
886 tilt_sw_output_types = Tuple.T(3, Int.T(), default=(0, 0, 0))
887 gravity_sw_output_types = Tuple.T(2, Int.T(), default=(0, 0))
889 indispl_filenames = Tuple.T(3, String.T(), default=psgrn_displ_names)
890 instress_filenames = Tuple.T(6, String.T(), default=psgrn_stress_names)
891 intilt_filenames = Tuple.T(3, String.T(), default=psgrn_tilt_names)
892 ingravity_filenames = Tuple.T(2, String.T(), default=psgrn_gravity_names)
894 outdispl_filenames = Tuple.T(3, String.T(), default=pscmp_displ_names)
895 outstress_filenames = Tuple.T(6, String.T(), default=pscmp_stress_names)
896 outtilt_filenames = Tuple.T(3, String.T(), default=pscmp_tilt_names)
897 outgravity_filenames = Tuple.T(2, String.T(), default=pscmp_gravity_names)
899 snapshot_basefilename = String.T(default='snapshot')
901 @classmethod
902 def example(cls):
903 conf = cls()
904 conf.psgrn_outdir = 'TEST_psgrn_functions/'
905 conf.pscmp_outdir = 'TEST_pscmp_output/'
906 conf.rectangular_source_patches = [PsCmpRectangularSource(
907 lat=10., lon=10., slip=2.,
908 width=5., length=10.,
909 strike=45, dip=30, rake=-90)]
910 conf.observation = PsCmpArray(
911 slat=9.5, elat=10.5, n_steps_lat=150,
912 slon=9.5, elon=10.5, n_steps_lon=150)
913 return conf
915 def get_output_filenames(self, rundir):
916 return [pjoin(rundir,
917 self.snapshot_basefilename + '_' + str(i + 1) + '.txt')
918 for i in range(len(self.snapshots.times))]
920 def string_for_config(self):
921 d = self.__dict__.copy()
923 patches_str, n_patches = distributed_fault_patches_to_config(
924 self.rectangular_source_patches)
926 d['patches_str'] = patches_str
927 d['n_patches'] = n_patches
929 str_npoints, str_observation = self.observation.string_for_config()
930 d['str_npoints'] = str_npoints
931 d['str_observation'] = str_observation
932 d['sw_observation_type'] = self.observation.sw
934 if self.snapshots.times:
935 str_times_dummy = []
936 for i, time in enumerate(self.snapshots.times):
937 str_times_dummy.append("%f '%s_%i.txt'\n" % (
938 time, self.snapshot_basefilename, i + 1))
940 str_times_dummy.append('#')
941 d['str_times_snapshots'] = ''.join(str_times_dummy)
942 d['n_snapshots'] = len(str_times_dummy) - 1
943 else:
944 d['str_times_snapshots'] = '# '
945 d['n_snapshots'] = 0
947 if self.sw_los_displacement == 1:
948 d['str_los_vector'] = str_float_vals(self.los_vector)
949 else:
950 d['str_los_vector'] = ''
952 if self.sw_coulomb_stress == 1:
953 d['str_coulomb_master_field'] = \
954 self.coulomb_master_field.string_for_config()
955 else:
956 d['str_coulomb_master_field'] = ''
958 d['str_psgrn_outdir'] = "'%s'" % self.psgrn_outdir
959 d['str_pscmp_outdir'] = "'%s'" % './'
961 d['str_indispl_filenames'] = str_str_vals(self.indispl_filenames)
962 d['str_instress_filenames'] = str_str_vals(self.instress_filenames)
963 d['str_intilt_filenames'] = str_str_vals(self.intilt_filenames)
964 d['str_ingravity_filenames'] = str_str_vals(self.ingravity_filenames)
966 d['str_outdispl_filenames'] = str_str_vals(self.outdispl_filenames)
967 d['str_outstress_filenames'] = str_str_vals(self.outstress_filenames)
968 d['str_outtilt_filenames'] = str_str_vals(self.outtilt_filenames)
969 d['str_outgravity_filenames'] = str_str_vals(self.outgravity_filenames)
971 d['str_displ_sw_output_types'] = str_int_vals(
972 self.displ_sw_output_types)
973 d['str_stress_sw_output_types'] = str_int_vals(
974 self.stress_sw_output_types)
975 d['str_tilt_sw_output_types'] = str_int_vals(
976 self.tilt_sw_output_types)
977 d['str_gravity_sw_output_types'] = str_int_vals(
978 self.gravity_sw_output_types)
980 template = '''# autogenerated PSCMP input by pscmp.py
981#===============================================================================
982# This is input file of FORTRAN77 program "pscmp08" for modeling post-seismic
983# deformation induced by earthquakes in multi-layered viscoelastic media using
984# the Green's function approach. The earthquke source is represented by an
985# arbitrary number of rectangular dislocation planes. For more details, please
986# read the accompanying READ.ME file.
987#
988# written by Rongjiang Wang
989# GeoForschungsZentrum Potsdam
990# e-mail: wang@gfz-potsdam.de
991# phone +49 331 2881209
992# fax +49 331 2881204
993#
994# Last modified: Potsdam, July, 2008
995#
996# References:
997#
998# (1) Wang, R., F. Lorenzo-Martin and F. Roth (2003), Computation of deformation
999# induced by earthquakes in a multi-layered elastic crust - FORTRAN programs
1000# EDGRN/EDCMP, Computer and Geosciences, 29(2), 195-207.
1001# (2) Wang, R., F. Lorenzo-Martin and F. Roth (2006), PSGRN/PSCMP - a new code for
1002# calculating co- and post-seismic deformation, geoid and gravity changes
1003# based on the viscoelastic-gravitational dislocation theory, Computers and
1004# Geosciences, 32, 527-541. DOI:10.1016/j.cageo.2005.08.006.
1005# (3) Wang, R. (2005), The dislocation theory: a consistent way for including the
1006# gravity effect in (visco)elastic plane-earth models, Geophysical Journal
1007# International, 161, 191-196.
1008#
1009#################################################################
1010## ##
1011## Green's functions should have been prepared with the ##
1012## program "psgrn08" before the program "pscmp08" is started. ##
1013## ##
1014## For local Cartesian coordinate system, the Aki's convention ##
1015## is used, that is, x is northward, y is eastward, and z is ##
1016## downward. ##
1017## ##
1018## If not specified otherwise, SI Unit System is used overall! ##
1019## ##
1020#################################################################
1021#===============================================================================
1022# OBSERVATION ARRAY
1023# =================
1024# 1. selection for irregular observation positions (= 0) or a 1D observation
1025# profile (= 1) or a rectangular 2D observation array (= 2): iposrec
1026#
1027# IF (iposrec = 0 for irregular observation positions) THEN
1028#
1029# 2. number of positions: nrec
1030#
1031# 3. coordinates of the observations: (lat(i),lon(i)), i=1,nrec
1032#
1033# ELSE IF (iposrec = 1 for regular 1D observation array) THEN
1034#
1035# 2. number of position samples of the profile: nrec
1036#
1037# 3. the start and end positions: (lat1,lon1), (lat2,lon2)
1038#
1039# ELSE IF (iposrec = 2 for rectanglular 2D observation array) THEN
1040#
1041# 2. number of x samples, start and end values: nxrec, xrec1, xrec2
1042#
1043# 3. number of y samples, start and end values: nyrec, yrec1, yrec2
1044#
1045# sequence of the positions in output data: lat(1),lon(1); ...; lat(nx),lon(1);
1046# lat(1),lon(2); ...; lat(nx),lon(2); ...; lat(1),lon(ny); ...; lat(nx),lon(ny).
1047#
1048# Note that the total number of observation positions (nrec or nxrec*nyrec)
1049# should be <= NRECMAX (see pecglob.h)!
1050#===============================================================================
1051 %(sw_observation_type)i
1052%(str_npoints)s
1053%(str_observation)s
1054#===============================================================================
1055# OUTPUTS
1056# =======
1057#
1058# 1. select output for los displacement (only for snapshots, see below), x, y,
1059# and z-cosines to the INSAR orbit: insar (1/0 = yes/no), xlos, ylos, zlos
1060#
1061# if this option is selected (insar = 1), the snapshots will include additional
1062# data:
1063# LOS_Dsp = los displacement to the given satellite orbit.
1064#
1065# 2. select output for Coulomb stress changes (only for snapshots, see below):
1066# icmb (1/0 = yes/no), friction, Skempton ratio, strike, dip, and rake angles
1067# [deg] describing the uniform regional master fault mechanism, the uniform
1068# regional principal stresses: sigma1, sigma2 and sigma3 [Pa] in arbitrary
1069# order (the orietation of the pre-stress field will be derived by assuming
1070# that the master fault is optimally oriented according to Coulomb failure
1071# criterion)
1072#
1073# if this option is selected (icmb = 1), the snapshots will include additional
1074# data:
1075# CMB_Fix, Sig_Fix = Coulomb and normal stress changes on master fault;
1076# CMB_Op1/2, Sig_Op1/2 = Coulomb and normal stress changes on the two optimally
1077# oriented faults;
1078# Str_Op1/2, Dip_Op1/2, Slp_Op1/2 = strike, dip and rake angles of the two
1079# optimally oriented faults.
1080#
1081# Note: the 1. optimally orieted fault is the one closest to the master fault.
1082#
1083# 3. output directory in char format: outdir
1084#
1085# 4. select outputs for displacement components (1/0 = yes/no): itout(i), i=1-3
1086#
1087# 5. the file names in char format for the x, y, and z components:
1088# toutfile(i), i=1-3
1089#
1090# 6. select outputs for stress components (1/0 = yes/no): itout(i), i=4-9
1091#
1092# 7. the file names in char format for the xx, yy, zz, xy, yz, and zx components:
1093# toutfile(i), i=4-9
1094#
1095# 8. select outputs for vertical NS and EW tilt components, block rotation, geoid
1096# and gravity changes (1/0 = yes/no): itout(i), i=10-14
1097#
1098# 9. the file names in char format for the NS tilt (positive if borehole top
1099# tilts to north), EW tilt (positive if borehole top tilts to east), block
1100# rotation (clockwise positive), geoid and gravity changes: toutfile(i), i=10-14
1101#
1102# Note that all above outputs are time series with the time window as same
1103# as used for the Green's functions
1104#
1105#10. number of scenario outputs ("snapshots": spatial distribution of all above
1106# observables at given time points; <= NSCENMAX (see pscglob.h): nsc
1107#
1108#11. the time [day], and file name (in char format) for the 1. snapshot;
1109#12. the time [day], and file name (in char format) for the 2. snapshot;
1110#13. ...
1111#
1112# Note that all file or directory names should not be longer than 80
1113# characters. Directories must be ended by / (unix) or \ (dos)!
1114#===============================================================================
1115 %(sw_los_displacement)i %(str_los_vector)s
1116 %(sw_coulomb_stress)i %(str_coulomb_master_field)s
1117 %(str_pscmp_outdir)s
1118 %(str_displ_sw_output_types)s
1119 %(str_outdispl_filenames)s
1120 %(str_stress_sw_output_types)s
1121 %(str_outstress_filenames)s
1122 %(str_tilt_sw_output_types)s %(str_gravity_sw_output_types)s
1123 %(str_outtilt_filenames)s %(str_outgravity_filenames)s
1124 %(n_snapshots)i
1125%(str_times_snapshots)s
1126#===============================================================================
1127#
1128# GREEN'S FUNCTION DATABASE
1129# =========================
1130# 1. directory where the Green's functions are stored: grndir
1131#
1132# 2. file names (without extensions!) for the 13 Green's functions:
1133# 3 displacement komponents (uz, ur, ut): green(i), i=1-3
1134# 6 stress components (szz, srr, stt, szr, srt, stz): green(i), i=4-9
1135# radial and tangential components measured by a borehole tiltmeter,
1136# rigid rotation around z-axis, geoid and gravity changes (tr, tt, rot, gd, gr):
1137# green(i), i=10-14
1138#
1139# Note that all file or directory names should not be longer than 80
1140# characters. Directories must be ended by / (unix) or \ (dos)! The
1141# extensions of the file names will be automatically considered. They
1142# are ".ep", ".ss", ".ds" and ".cl" denoting the explosion (inflation)
1143# strike-slip, the dip-slip and the compensated linear vector dipole
1144# sources, respectively.
1145#
1146#===============================================================================
1147 %(str_psgrn_outdir)s
1148 %(str_indispl_filenames)s
1149 %(str_instress_filenames)s
1150 %(str_intilt_filenames)s %(str_ingravity_filenames)s
1151#===============================================================================
1152# RECTANGULAR SUBFAULTS
1153# =====================
1154# 1. number of subfaults (<= NSMAX in pscglob.h): ns
1155#
1156# 2. parameters for the 1. rectangular subfault: geographic coordinates
1157# (O_lat, O_lon) [deg] and O_depth [km] of the local reference point on
1158# the present fault plane, length (along strike) [km] and width (along down
1159# dip) [km], strike [deg], dip [deg], number of equi-size fault patches along
1160# the strike (np_st) and along the dip (np_di) (total number of fault patches
1161# = np_st x np_di), and the start time of the rupture; the following data
1162# lines describe the slip distribution on the present sub-fault:
1163#
1164# pos_s[km] pos_d[km] slip_strike[m] slip_downdip[m] opening[m]
1165#
1166# where (pos_s,pos_d) defines the position of the center of each patch in
1167# the local coordinate system with the origin at the reference point:
1168# pos_s = distance along the length (positive in the strike direction)
1169# pos_d = distance along the width (positive in the down-dip direction)
1170#
1171#
1172# 3. ... for the 2. subfault ...
1173# ...
1174# N
1175# /
1176# /| strike
1177# +------------------------
1178# |\ p . \ W
1179# :-\ i . \ i
1180# | \ l . \ d
1181# :90 \ S . \ t
1182# |-dip\ . \ h
1183# : \. | rake \
1184# Z -------------------------
1185# L e n g t h
1186#
1187# Simulation of a Mogi source:
1188# (1) Calculate deformation caused by three small openning plates (each
1189# causes a third part of the volume of the point inflation) located
1190# at the same depth as the Mogi source but oriented orthogonal to
1191# each other.
1192# (2) Multiply the results by 3(1-nu)/(1+nu), where nu is the Poisson
1193# ratio at the source depth.
1194# The multiplication factor is the ratio of the seismic moment (energy) of
1195# the Mogi source to that of the plate openning with the same volume change.
1196#===============================================================================
1197# n_faults
1198#-------------------------------------------------------------------------------
1199 %(n_patches)i
1200#-------------------------------------------------------------------------------
1201# n O_lat O_lon O_depth length width strike dip np_st np_di start_time
1202# [-] [deg] [deg] [km] [km] [km] [deg] [deg] [-] [-] [day]
1203# pos_s pos_d slp_stk slp_ddip open
1204# [km] [km] [m] [m] [m]
1205#-------------------------------------------------------------------------------
1206%(patches_str)s
1207#================================end of input===================================
1208''' # noqa
1209 return (template % d).encode('ascii')
1212class PsGrnPsCmpConfig(Object):
1213 '''
1214 Combined config Object of PsGrn and PsCmp.
1216 Note: attributes in this class use non-standard units [days] to be
1217 consistent with PSCMP text file input. Please read the documentation
1218 carefully.
1219 '''
1220 tmin_days = Float.T(
1221 default=0.,
1222 help='Min. time in days')
1223 tmax_days = Float.T(
1224 default=1.,
1225 help='Max. time in days')
1226 gf_outdir = String.T(default='psgrn_functions')
1228 psgrn_config = PsGrnConfig.T(default=PsGrnConfig.D())
1229 pscmp_config = PsCmpConfig.T(default=PsCmpConfig.D())
1232class PsCmpError(gf.store.StoreError):
1233 pass
1236class Interrupted(gf.store.StoreError):
1237 def __str__(self):
1238 return 'Interrupted.'
1241class PsCmpRunner(object):
1242 '''
1243 Wrapper object to execute the program fomosto_pscmp with the specified
1244 configuration.
1246 :param tmp: string, path to the temporary directy where calculation
1247 results are stored
1248 :param keep_tmp: boolean, if True the result directory is kept
1249 '''
1250 def __init__(self, tmp=None, keep_tmp=False):
1251 if tmp is not None:
1252 tmp = os.path.abspath(tmp)
1253 self.tempdir = mkdtemp(prefix='pscmprun-', dir=tmp)
1254 self.keep_tmp = keep_tmp
1255 self.config = None
1257 def run(self, config):
1258 '''
1259 Run the program!
1261 :param config: :py:class:`PsCmpConfigFull`
1262 '''
1263 self.config = config
1265 input_fn = pjoin(self.tempdir, 'input')
1267 with open(input_fn, 'wb') as f:
1268 input_str = config.string_for_config()
1270 logger.debug('===== begin pscmp input =====\n'
1271 '%s===== end pscmp input =====' % input_str.decode())
1273 f.write(input_str)
1275 program = program_bins['pscmp.%s' % config.version]
1277 old_wd = os.getcwd()
1278 os.chdir(self.tempdir)
1280 interrupted = []
1282 def signal_handler(signum, frame):
1283 os.kill(proc.pid, signal.SIGTERM)
1284 interrupted.append(True)
1286 original = signal.signal(signal.SIGINT, signal_handler)
1287 try:
1288 try:
1289 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE,
1290 close_fds=True)
1292 except OSError as err:
1293 os.chdir(old_wd)
1294 logger.error('OS error: {0}'.format(err))
1295 raise PsCmpError(
1296 '''could not start pscmp executable: "%s"
1297Available fomosto backends and download links to the modelling codes are listed
1298on
1300 https://pyrocko.org/docs/current/apps/fomosto/backends.html
1302''' % program)
1304 (output_str, error_str) = proc.communicate(b'input\n')
1306 finally:
1307 signal.signal(signal.SIGINT, original)
1309 if interrupted:
1310 raise KeyboardInterrupt()
1312 logger.debug('===== begin pscmp output =====\n'
1313 '%s===== end pscmp output =====' % output_str.decode())
1315 errmsg = []
1316 if proc.returncode != 0:
1317 errmsg.append(
1318 'pscmp had a non-zero exit state: %i' % proc.returncode)
1320 if error_str:
1321 errmsg.append('pscmp emitted something via stderr')
1323 if output_str.lower().find(b'error') != -1:
1324 errmsg.append("the string 'error' appeared in pscmp output")
1326 if errmsg:
1327 self.keep_tmp = True
1329 os.chdir(old_wd)
1330 raise PsCmpError('''
1331===== begin pscmp input =====
1332{pscmp_input}===== end pscmp input =====
1333===== begin pscmp output =====
1334{pscmp_output}===== end pscmp output =====
1335===== begin pscmp error =====
1336{pscmp_error}===== end pscmp error =====
1337{error_messages}
1338pscmp has been invoked as "{call}"
1339in the directory {dir}'''.format(
1340 pscmp_input=input_str,
1341 pscmp_output=output_str,
1342 pscmp_error=error_str,
1343 error_messages='\n'.join(errmsg),
1344 call=program,
1345 dir=self.tempdir)
1346 .lstrip())
1348 self.pscmp_output = output_str
1349 self.pscmp_error = error_str
1351 os.chdir(old_wd)
1353 def get_results(self, component='displ'):
1354 '''
1355 Get the resulting components from the stored directory.
1357 .. note::
1359 The z-component is downward positive.
1361 :param component:
1362 The component to retrieve from the result directory. Choices:
1363 ``"displ"`` -- displacement, n x 3 array,
1364 ``"stress"`` -- stresses n x 6 array,
1365 ``"tilt'`` -- tilts n x 3 array,
1366 ``"gravity'`` -- gravity n x 2 array,
1367 ``"all"`` -- all the above together.
1368 :type param:
1369 str
1370 '''
1371 assert self.config.snapshots is not None
1372 fns = self.config.get_output_filenames(self.tempdir)
1374 output = []
1375 for fn in fns:
1376 if not os.path.exists(fn):
1377 continue
1379 data = num.loadtxt(fn, skiprows=1, dtype=float)
1381 try:
1382 _, idxs = pscmp_component_mapping[component]
1383 except KeyError:
1384 raise Exception('component %s not supported! Either: %s' % (
1385 component, ', '.join(
1386 '"%s"' % k for k in pscmp_component_mapping.keys())))
1388 istart, iend = idxs
1390 output.append(data[:, istart:iend])
1392 return output
1394 def get_traces(self, component='displ'):
1395 '''
1396 Load snapshot data array and return specified components.
1397 Transform array component and receiver wise to list of
1398 :py:class:`pyrocko.trace.Trace`
1399 '''
1401 distances = self.config.observation.distances
1403 output_list = self.get_results(component=component)
1404 deltat = self.config.snapshots.deltat
1406 nrec, ncomp = output_list[0].shape
1408 outarray = num.dstack([num.zeros([nrec, ncomp])] + output_list)
1410 comps, _ = pscmp_component_mapping[component]
1412 outtraces = []
1413 for row in range(nrec):
1414 for col, comp in enumerate(comps):
1415 tr = trace.Trace(
1416 '', '%04i' % row, '', comp,
1417 tmin=0.0, deltat=deltat, ydata=outarray[row, col, :],
1418 meta=dict(distance=distances[row]))
1419 outtraces.append(tr)
1421 return outtraces
1423 def __del__(self):
1424 if self.tempdir:
1425 if not self.keep_tmp:
1426 shutil.rmtree(self.tempdir)
1427 self.tempdir = None
1428 else:
1429 logger.warning(
1430 'not removing temporary directory: %s' % self.tempdir)
1433class PsGrnCmpGFBuilder(gf.builder.Builder):
1434 nsteps = 2
1436 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
1437 force=False):
1439 self.store = gf.store.Store(store_dir, 'w')
1441 storeconf = self.store.config
1443 dummy_lat = 10.0
1444 dummy_lon = 10.0
1446 depths = storeconf.coords[0]
1447 lats = num.ones_like(depths) * dummy_lat
1448 lons = num.ones_like(depths) * dummy_lon
1449 points = num.vstack([lats, lons, depths]).T
1451 self.shear_moduli = storeconf.get_shear_moduli(
1452 lat=dummy_lat,
1453 lon=dummy_lon,
1454 points=points,
1455 interpolation='multilinear')
1457 self.lambda_moduli = storeconf.get_lambda_moduli(
1458 lat=dummy_lat,
1459 lon=dummy_lon,
1460 points=points,
1461 interpolation='multilinear')
1463 if step == 0:
1464 block_size = (1, storeconf.nsource_depths, storeconf.ndistances)
1465 else:
1466 if block_size is None:
1467 block_size = (1, 1, storeconf.ndistances)
1469 if len(storeconf.ns) == 2:
1470 block_size = block_size[1:]
1472 gf.builder.Builder.__init__(
1473 self, storeconf, step, block_size=block_size, force=force)
1475 baseconf = self.store.get_extra('psgrn_pscmp')
1477 cg = PsGrnConfigFull(**baseconf.psgrn_config.items())
1478 if cg.n_time_samples is None or cg.max_time is None:
1479 deltat_days = 1. / storeconf.sample_rate / day
1481 cg.n_time_samples = int(baseconf.tmax_days // deltat_days)
1482 cg.max_time = baseconf.tmax_days
1483 else:
1484 warnings.warn(
1485 'PsGrnConfig is defining n_times_samples and max_time,'
1486 ' this is replaced by PsGrnPsCmpConfig tmin and tmax.',
1487 FutureWarning)
1489 cc = PsCmpConfigFull(**baseconf.pscmp_config.items())
1490 if cc.snapshots is None:
1491 deltat_days = 1. / storeconf.sample_rate / day
1493 cc.snapshots = PsCmpSnapshots(
1494 tmin=baseconf.tmin_days,
1495 tmax=baseconf.tmax_days,
1496 deltatdays=deltat_days)
1497 else:
1498 warnings.warn(
1499 'PsCmpConfig is defining snapshots,'
1500 ' this is replaced by PsGrnPsCmpConfig tmin and tmax.',
1501 FutureWarning)
1503 cg.earthmodel_1d = storeconf.earthmodel_1d
1505 gf_outpath = os.path.join(store_dir, baseconf.gf_outdir)
1507 cg.psgrn_outdir = gf_outpath + '/'
1508 cc.psgrn_outdir = gf_outpath + '/'
1510 util.ensuredir(gf_outpath)
1512 self.cg = cg
1513 self.cc = cc
1515 def cleanup(self):
1516 self.store.close()
1518 def work_block(self, iblock):
1520 if len(self.store.config.ns) == 2:
1521 (sz, firstx), (sz, lastx), (ns, nx) = \
1522 self.get_block_extents(iblock)
1523 mu = self.shear_moduli[iblock]
1524 lame_lambda = self.lambda_moduli[iblock]
1526 rz = self.store.config.receiver_depth
1527 else:
1528 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
1529 self.get_block_extents(iblock)
1531 fc = self.store.config
1532 cc = copy.deepcopy(self.cc)
1533 cg = copy.deepcopy(self.cg)
1535 dx = fc.distance_delta
1537 logger.info(
1538 'Starting step %i / %i, block %i / %i' %
1539 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
1541 if self.step == 0:
1542 if cg.gf_depth_spacing is None:
1543 gf_depth_spacing = fc.source_depth_delta
1544 else:
1545 gf_depth_spacing = cg.gf_depth_spacing * km
1547 n_steps_depth = int((fc.source_depth_max - fc.source_depth_min) /
1548 gf_depth_spacing) + 1
1550 if cg.gf_distance_spacing is None:
1551 gf_distance_spacing = fc.distance_delta
1552 else:
1553 gf_distance_spacing = cg.gf_distance_spacing * km
1554 n_steps_distance = int((fc.distance_max - fc.distance_min) /
1555 gf_distance_spacing) + 1
1557 cg.depth_grid = PsGrnSpatialSampling(
1558 n_steps=n_steps_depth,
1559 start_distance=fc.source_depth_min / km,
1560 end_distance=fc.source_depth_max / km)
1562 cg.distance_grid = PsGrnSpatialSampling(
1563 n_steps=n_steps_distance,
1564 start_distance=fc.distance_min / km,
1565 end_distance=fc.distance_max / km)
1567 runner = PsGrnRunner(outdir=self.cg.psgrn_outdir)
1568 runner.run(cg, force=self.force)
1570 else:
1571 distances = num.linspace(
1572 firstx, firstx + (nx - 1) * dx, nx).tolist()
1574 # fomosto sample rate in s, pscmp takes days
1575 deltatdays = 1. / (fc.sample_rate * 24. * 3600.)
1576 cc.snapshots = PsCmpSnapshots(
1577 tmin=0., tmax=cg.max_time, deltatdays=deltatdays)
1578 cc.observation = PsCmpProfile(
1579 slat=0. - 0.001 * cake.m2d,
1580 slon=0.0,
1581 elat=0. + distances[-1] * cake.m2d,
1582 elon=0.0,
1583 n_steps=len(distances),
1584 distances=distances)
1586 runner = PsCmpRunner(keep_tmp=False)
1588 mtsize = float(dx * cc.rectangular_fault_size_factor)
1590 Ai = 1. / num.power(mtsize, 2)
1591 nullf, sf = get_iso_scaling_factors(mu=mu, lame_lambda=lame_lambda)
1593 mui = 1. / mu
1595 gfmapping = [
1596 (('nn',),
1597 {'un': (0, Ai * sf), 'ud': (5, Ai * sf)}),
1598 (('ne',),
1599 {'ue': (3, 1 * Ai * mui)}),
1600 (('nd',),
1601 {'un': (1, 1 * Ai * mui), 'ud': (6, 1 * Ai * mui)}),
1602 (('ed',),
1603 {'ue': (4, 1 * Ai * mui)}),
1604 (('dd',),
1605 {'un': (2, Ai * sf), 'ud': (7, Ai * sf)}),
1606 (('ee',),
1607 {'un': (8, Ai * sf), 'ud': (9, Ai * sf)}),
1608 ]
1610 for mt, gfmap in gfmapping:
1611 cc.rectangular_source_patches = []
1612 for idx in mt:
1613 pmt = PsCmpMomentTensor(
1614 lat=0. + 0.001 * dx * cake.m2d,
1615 lon=0.0,
1616 depth=float(sz),
1617 width=mtsize,
1618 length=mtsize,
1619 idx=idx)
1621 cc.rectangular_source_patches.extend(
1622 pmt.to_rfs(nullf))
1624 runner.run(cc)
1626 rawtraces = runner.get_traces()
1628 interrupted = []
1630 def signal_handler(signum, frame):
1631 interrupted.append(True)
1633 original = signal.signal(signal.SIGINT, signal_handler)
1634 self.store.lock()
1635 duplicate_inserts = 0
1637 try:
1638 for itr, tr in enumerate(rawtraces):
1639 if tr.channel in gfmap:
1640 x = tr.meta['distance']
1641 ig, factor = gfmap[tr.channel]
1643 if len(self.store.config.ns) == 2:
1644 args = (sz, x, ig)
1645 else:
1646 args = (rz, sz, x, ig)
1648 gf_tr = gf.store.GFTrace.from_trace(tr)
1650 gf_tr.data *= factor
1652 try:
1653 self.store.put(args, gf_tr)
1654 except gf.store.DuplicateInsert:
1655 duplicate_inserts += 1
1657 finally:
1658 if duplicate_inserts:
1659 logger.warning(
1660 '%i insertions skipped (duplicates)'
1661 % duplicate_inserts)
1663 self.store.unlock()
1664 signal.signal(signal.SIGINT, original)
1666 if interrupted:
1667 raise KeyboardInterrupt()
1669 logger.info(
1670 'Done with step %i / %i, block %i / %i' % (
1671 self.step + 1, self.nsteps, iblock + 1, self.nblocks))
1674def init(store_dir, variant):
1675 if variant is None:
1676 variant = '2008a'
1678 if ('pscmp.' + variant) not in program_bins:
1679 raise gf.store.StoreError('unsupported pscmp variant: %s' % variant)
1681 if ('psgrn.' + variant) not in program_bins:
1682 raise gf.store.StoreError('unsupported psgrn variant: %s' % variant)
1684 c = PsGrnPsCmpConfig()
1686 store_id = os.path.basename(os.path.realpath(store_dir))
1688 # Initialising a viscous mantle
1689 cake_mod = cake.load_model(fn=None, crust2_profile=(54., 23.))
1690 mantle = cake_mod.material(z=45*km)
1691 mantle.burger_eta1 = 5e17
1692 mantle.burger_eta2 = 1e19
1693 mantle.burger_alpha = 1.
1695 config = gf.meta.ConfigTypeA(
1696 id=store_id,
1697 ncomponents=10,
1698 sample_rate=1. / (3600. * 24.),
1699 receiver_depth=0. * km,
1700 source_depth_min=0. * km,
1701 source_depth_max=15. * km,
1702 source_depth_delta=.5 * km,
1703 distance_min=0. * km,
1704 distance_max=50. * km,
1705 distance_delta=1. * km,
1706 earthmodel_1d=cake_mod,
1707 modelling_code_id='psgrn_pscmp.%s' % variant,
1708 tabulated_phases=[]) # dummy list
1710 c.validate()
1711 config.validate()
1712 return gf.store.Store.create_editables(
1713 store_dir,
1714 config=config,
1715 extra={'psgrn_pscmp': c})
1718def build(store_dir,
1719 force=False,
1720 nworkers=None,
1721 continue_=False,
1722 step=None,
1723 iblock=None):
1725 return PsGrnCmpGFBuilder.build(
1726 store_dir, force=force, nworkers=nworkers, continue_=continue_,
1727 step=step, iblock=iblock)