1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import numpy as num
8import logging
9import os
10import shutil
11import glob
12import copy
13import signal
14import math
15import time
17from tempfile import mkdtemp
18from subprocess import Popen, PIPE
19from os.path import join as pjoin
21from pyrocko import trace, util, cake, gf
22from pyrocko.guts import Object, Float, String, Bool, Tuple, Int, List
23from pyrocko.moment_tensor import MomentTensor, symmat6
25guts_prefix = 'pf'
27logger = logging.getLogger('pyrocko.fomosto.qssp')
29# how to call the programs
30program_bins = {
31 'qssp.2010beta': 'fomosto_qssp2010beta',
32 'qssp.2010': 'fomosto_qssp2010',
33 'qssp.2017': 'fomosto_qssp2017',
34 'qssp.2020': 'fomosto_qssp2020',
35 'qssp.ppeg2017': 'fomosto_qsspppeg2017',
36}
39def have_backend(version=None):
40 avail = set()
41 for version, exe in program_bins.items():
42 try:
43 p = Popen([exe], stdout=PIPE, stderr=PIPE, stdin=PIPE)
44 (stdout, stderr) = p.communicate()
45 avail.add(version)
47 except OSError:
48 pass
50 return avail
53qssp_components = {
54 1: 'ae an az gr sd te tn ue un uz ve vn vz'.split(),
55 2: 'ar at ap gr sd tt tp ur ut up vr vt vp'.split(),
56 3: '_disp_e _disp_n _disp_z'.split(),
57 4: '_gravitation_e _gravitation_n _gravitation_z '
58 '_acce_e _acce_n _acce_z'.split(),
59 5: '_rota_e _rota_n _rota_z'.split()
60}
63def str_float_vals(vals):
64 return ' '.join(['%.6e' % val for val in vals])
67def cake_model_to_config(mod):
68 k = 1000.
69 srows = []
70 for i, row in enumerate(mod.to_scanlines()):
71 depth, vp, vs, rho, qp, qs = row
72 row = [depth/k, vp/k, vs/k, rho/k, qp, qs]
73 srows.append('%i %s' % (i+1, str_float_vals(row)))
75 return '\n'.join(srows), len(srows)
78class QSSPSource(Object):
79 lat = Float.T(default=0.0)
80 lon = Float.T(default=0.0)
81 depth = Float.T(default=10.0)
82 torigin = Float.T(default=0.0)
83 trise = Float.T(default=1.0)
85 def string_for_config(self):
86 return '%(lat)15e %(lon)15e %(depth)15e %(torigin)15e %(trise)15e' \
87 % self.__dict__
90class QSSPSourceMT(QSSPSource):
91 munit = Float.T(default=1.0)
92 mrr = Float.T(default=1.0)
93 mtt = Float.T(default=1.0)
94 mpp = Float.T(default=1.0)
95 mrt = Float.T(default=0.0)
96 mrp = Float.T(default=0.0)
97 mtp = Float.T(default=0.0)
99 def string_for_config(self):
100 return '%(munit)15e %(mrr)15e %(mtt)15e %(mpp)15e ' \
101 '%(mrt)15e %(mrp)15e %(mtp)15e ' \
102 % self.__dict__ + QSSPSource.string_for_config(self)
105class QSSPSourceDC(QSSPSource):
106 moment = Float.T(default=1.0e9)
107 strike = Float.T(default=0.0)
108 dip = Float.T(default=90.0)
109 rake = Float.T(default=0.0)
111 def string_for_config(self):
112 return '%(moment)15e %(strike)15e %(dip)15e %(rake)15e ' \
113 % self.__dict__ + QSSPSource.string_for_config(self)
116class QSSPReceiver(Object):
117 lat = Float.T(default=10.0)
118 lon = Float.T(default=0.0)
119 name = String.T(default='')
120 tstart = Float.T(default=0.0)
121 distance = Float.T(default=0.0)
123 def string_for_config(self):
124 return "%(lat)15e %(lon)15e '%(name)s' %(tstart)e" % self.__dict__
127class QSSPGreen(Object):
128 depth = Float.T(default=10.0)
129 filename = String.T(default='GF_10km')
130 calculate = Bool.T(default=True)
132 def string_for_config(self):
133 return "%(depth)15e '%(filename)s' %(calculate)i" % self.__dict__
136class QSSPConfig(Object):
137 qssp_version = String.T(default='2010beta')
138 time_region = Tuple.T(2, gf.Timing.T(), default=(
139 gf.Timing('-10'), gf.Timing('+890')))
141 frequency_max = Float.T(optional=True)
142 slowness_max = Float.T(default=0.4)
143 antialiasing_factor = Float.T(default=0.1)
145 # only available in 2017:
146 switch_turning_point_filter = Int.T(default=0)
147 max_pene_d1 = Float.T(default=2891.5)
148 max_pene_d2 = Float.T(default=6371.0)
149 earth_radius = Float.T(default=6371.0)
150 switch_free_surf_reflection = Int.T(default=1)
152 lowpass_order = Int.T(default=0, optional=True)
153 lowpass_corner = Float.T(default=1.0, optional=True)
155 bandpass_order = Int.T(default=0, optional=True)
156 bandpass_corner_low = Float.T(default=1.0, optional=True)
157 bandpass_corner_high = Float.T(default=1.0, optional=True)
159 output_slowness_min = Float.T(default=0.0, optional=True)
160 output_slowness_max = Float.T(optional=True)
162 spheroidal_modes = Bool.T(default=True)
163 toroidal_modes = Bool.T(default=True)
165 # only available in 2010beta:
166 cutoff_harmonic_degree_sd = Int.T(optional=True, default=0)
168 cutoff_harmonic_degree_min = Int.T(default=0)
169 cutoff_harmonic_degree_max = Int.T(default=25000)
171 crit_frequency_sge = Float.T(default=0.0)
172 crit_harmonic_degree_sge = Int.T(default=0)
174 include_physical_dispersion = Bool.T(default=False)
176 source_patch_radius = Float.T(default=0.0)
178 cut = Tuple.T(2, gf.Timing.T(), optional=True)
180 fade = Tuple.T(4, gf.Timing.T(), optional=True)
181 relevel_with_fade_in = Bool.T(default=False)
182 nonzero_fade_in = Bool.T(default=False)
183 nonzero_fade_out = Bool.T(default=False)
185 def items(self):
186 return dict(self.T.inamevals(self))
189class QSSPConfigFull(QSSPConfig):
190 time_window = Float.T(default=900.0)
192 receiver_depth = Float.T(default=0.0)
193 sampling_interval = Float.T(default=5.0)
195 stored_quantity = String.T(default='displacement')
196 output_filename = String.T(default='receivers')
197 output_format = Int.T(default=1)
198 output_time_window = Float.T(optional=True)
200 gf_directory = String.T(default='qssp_green')
201 greens_functions = List.T(QSSPGreen.T())
203 sources = List.T(QSSPSource.T())
204 receivers = List.T(QSSPReceiver.T())
206 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True)
208 @staticmethod
209 def example():
210 conf = QSSPConfigFull()
211 conf.sources.append(QSSPSourceMT())
212 lats = [20.]
213 conf.receivers.extend(QSSPReceiver(lat=lat) for lat in lats)
214 conf.greens_functions.append(QSSPGreen())
215 return conf
217 @property
218 def components(self):
219 if self.qssp_version in ('2017', '2020'):
220 if self.stored_quantity == "rotation":
221 fmt = 5
222 else:
223 fmt = 3
224 elif self.qssp_version == 'ppeg2017':
225 fmt = 4
226 else:
227 fmt = self.output_format
229 return qssp_components[fmt]
231 def get_output_filenames(self, rundir):
232 if self.qssp_version in ('2017', '2020', 'ppeg2017'):
233 return [
234 pjoin(rundir, self.output_filename + c + '.dat')
235 for c in self.components]
236 else:
237 return [
238 pjoin(rundir, self.output_filename + '.' + c)
239 for c in self.components]
241 def ensure_gf_directory(self):
242 util.ensuredir(self.gf_directory)
244 def string_for_config(self):
246 def aggregate(xx):
247 return len(xx), '\n'.join(x.string_for_config() for x in xx)
249 assert len(self.greens_functions) > 0
250 assert len(self.sources) > 0
251 assert len(self.receivers) > 0
253 d = self.__dict__.copy()
255 if self.output_time_window is None:
256 d['output_time_window'] = self.time_window
258 if self.output_slowness_max is None:
259 d['output_slowness_max'] = self.slowness_max
261 if self.frequency_max is None:
262 d['frequency_max'] = 0.5/self.sampling_interval
264 d['gf_directory'] = os.path.abspath(self.gf_directory) + '/'
266 d['n_receiver_lines'], d['receiver_lines'] = aggregate(self.receivers)
267 d['n_source_lines'], d['source_lines'] = aggregate(self.sources)
268 d['n_gf_lines'], d['gf_lines'] = aggregate(self.greens_functions)
269 model_str, nlines = cake_model_to_config(self.earthmodel_1d)
270 d['n_model_lines'] = nlines
271 d['model_lines'] = model_str
272 if self.stored_quantity == "rotation":
273 d['output_rotation'] = 1
274 d['output_displacement'] = 0
275 else:
276 d['output_displacement'] = 1
277 d['output_rotation'] = 0
279 if len(self.sources) == 0 or isinstance(self.sources[0], QSSPSourceMT):
280 d['point_source_type'] = 1
281 else:
282 d['point_source_type'] = 2
284 if self.qssp_version == '2010beta':
285 d['scutoff_doc'] = '''
286# (SH waves), and cutoff harmonic degree for static deformation
287'''.strip()
289 d['scutoff'] = '%i' % self.cutoff_harmonic_degree_sd
291 d['sfilter_doc'] = '''
292# 3. selection of order of Butterworth low-pass filter (if <= 0, then no
293# filtering), corner frequency (smaller than the cut-off frequency defined
294# above)
295'''.strip()
297 if self.bandpass_order != 0:
298 raise QSSPError(
299 'this version of qssp does not support bandpass '
300 'settings, use lowpass instead')
302 d['sfilter'] = '%i %f' % (
303 self.lowpass_order,
304 self.lowpass_corner)
306 elif self.qssp_version in ('2010', '2017', '2020', 'ppeg2017'):
307 d['scutoff_doc'] = '''
308# (SH waves), minimum and maximum cutoff harmonic degrees
309# Note: if the near-field static displacement is desired, the minimum
310# cutoff harmonic degree should not be smaller than, e.g., 2000.
311'''.strip()
313 d['scutoff'] = '%i %i' % (
314 self.cutoff_harmonic_degree_min,
315 self.cutoff_harmonic_degree_max)
317 d['sfilter_doc'] = '''
318# 3. selection of order of Butterworth bandpass filter (if <= 0, then no
319# filtering), lower and upper corner frequencies (smaller than the cut-off
320# frequency defined above)
321'''.strip()
323 if self.lowpass_order != 0:
324 raise QSSPError(
325 'this version of qssp does not support lowpass settings, '
326 'use bandpass instead')
328 d['sfilter'] = '%i %f %f' % (
329 self.bandpass_order,
330 self.bandpass_corner_low,
331 self.bandpass_corner_high)
332 if self.qssp_version in ('2017', '2020', 'ppeg2017'):
333 template = '''# autogenerated QSSP input by qssp.py
334#
335# This is the input file of FORTRAN77 program "qssp2017" for calculating
336# synthetic seismograms of a self-gravitating, spherically symmetric,
337# isotropic and viscoelastic earth.
338#
339# by
340# Rongjiang Wang <wang@gfz-potsdam.de>
341# Helmholtz-Centre Potsdam
342# GFZ German Reseach Centre for Geosciences
343# Telegrafenberg, D-14473 Potsdam, Germany
344#
345# Last modified: Potsdam, October 2017
346#
347# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
348# If not specified, SI Unit System is used overall!
349#
350# Coordinate systems:
351# spherical (r,t,p) with r = radial,
352# t = co-latitude,
353# p = east longitude.
354# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
355#
356# UNIFORM RECEIVER DEPTH
357# ======================
358# 1. uniform receiver depth [km]
359#-------------------------------------------------------------------------------------------
360 %(receiver_depth)e
361#-------------------------------------------------------------------------------------------
362#
363# SPACE-TIME SAMPLING PARAMETERS
364# =========================
365# 1. time window [sec], sampling interval [sec]
366# 2. max. frequency [Hz] of Green's functions
367# 3. max. slowness [s/km] of Green's functions
368# Note: if the near-field static displacement is desired, the maximum slowness should not
369# be smaller than the S wave slowness in the receiver layer
370# 4. anti-aliasing factor (> 0 & < 1), if it is <= 0 or >= 1/e (~ 0.4), then
371# default value of 1/e is used (e.g., 0.1 = alias phases will be suppressed
372# to 10%% of their original amplitude)
373# 5. switch (1/0 = yes/no) of turning-point filter, the range (d1, d2) of max. penetration
374# depth [km] (d1 is meaningless if it is smaller than the receiver/source depth, and
375# d2 is meaningless if it is equal to or larger than the earth radius)
376#
377# Note: The turning-point filter (Line 5) works only for the extended QSSP code (e.g.,
378# qssp2016). if this filter is selected, all phases with the turning point
379# shallower than d1 or deeper than d2 will be filtered.
380#
381# 6. Earth radius [km], switch of free-surface-reflection filter (1/0 = with/without free
382# surface reflection)
383#
384# Note: The free-surface-reflection filter (Line 6) works only for the extended QSSP
385# code (e.g., qssp2016). if this filter is selected, all phases with the turning
386# point shallower than d1 or deeper than d2 will be filtered.
387#-------------------------------------------------------------------------------------------
388 %(time_window)e %(sampling_interval)e
389 %(frequency_max)e
390 %(slowness_max)e
391 %(antialiasing_factor)e
392 %(switch_turning_point_filter)i %(max_pene_d1)e %(max_pene_d2)e
393 %(earth_radius)e %(switch_free_surf_reflection)i
394#-------------------------------------------------------------------------------------------
395#
396# SELF-GRAVITATING EFFECT
397# =======================
398# 1. the critical frequency [Hz] and the critical harmonic degree, below which
399# the self-gravitating effect should be included
400#-------------------------------------------------------------------------------------------
401 %(crit_frequency_sge)e %(crit_harmonic_degree_sge)i
402#-------------------------------------------------------------------------------------------
403#
404# WAVE TYPES
405# ==========
406# 1. selection (1/0 = yes/no) of speroidal modes (P-SV waves), selection of toroidal modes
407%(scutoff_doc)s
408#-------------------------------------------------------------------------------------------
409 %(spheroidal_modes)i %(toroidal_modes)i %(scutoff)s
410#-------------------------------------------------------------------------------------------
411# GREEN'S FUNCTION FILES
412# ======================
413# 1. number of discrete source depths, estimated radius of each source patch [km] and
414# directory for Green's functions
415# 2. list of the source depths [km], the respective file names of the Green's
416# functions (spectra) and the switch number (0/1) (0 = do not calculate
417# this Green's function because it exists already, 1 = calculate or update
418# this Green's function. Note: update is required if any of the above
419# parameters is changed)
420#-------------------------------------------------------------------------------------------
421 %(n_gf_lines)i %(source_patch_radius)e '%(gf_directory)s'
422 %(gf_lines)s
423#--------------------------------------------------------------------------------------------------------
424#
425# MULTI-EVENT SOURCE PARAMETERS
426# =============================
427# 1. number of discrete point sources and selection of the source data format
428# (1, 2 or 3)
429# 2. list of the multi-event sources
430#
431# Format 1 (full moment tensor):
432# Unit Mrr Mtt Mpp Mrt Mrp Mtp Lat Lon Depth T_origin T_rise
433# [Nm] [deg] [deg] [km] [sec] [sec]
434#
435# Format 2 (double couple):
436# Unit Strike Dip Rake Lat Lon Depth T_origin T_rise
437# [Nm] [deg] [deg] [deg] [deg] [deg] [km] [sec] [sec]
438#
439# Format 3 (single force):
440# Unit Feast Fnorth Fvertical Lat Lon Depth T_origin T_rise
441# [N] [deg] [deg] [km] [sec] [sec]
442#
443# Note: for each point source, the default moment (force) rate time function is used, defined by a
444# squared half-period (T_rise) sinusoid starting at T_origin.
445#-----------------------------------------------------------------------------------
446 %(n_source_lines)i %(point_source_type)i
447%(source_lines)s
448#--------------------------------------------------------------------------------------------------------
449#
450# RECEIVER PARAMETERS
451# ===================
452# 1. select output observables (1/0 = yes/no)
453# Note: the gravity change defined here is space based, i.e., the effect due to free-air
454# gradient and inertial are not included. the vertical component is positve upwards.
455# 2. output file name
456# 3. output time window [sec] (<= Green's function time window)
457%(sfilter_doc)s
458# 5. lower and upper slowness cut-off [s/km] (slowness band-pass filter)
459# 6. number of receiver
460# 7. list of the station parameters
461# Format:
462# Lat Lon Name Time_reduction
463# [deg] [deg] [sec]
464# (Note: Time_reduction = start time of the time window)
465#---------------------------------------------------------------------------------------------------------------
466# disp | velo | acce | strain | strain_rate | stress | stress_rate | rotation | rot_rate | gravitation | gravity
467#---------------------------------------------------------------------------------------------------------------
468# 1 1 1 1 1 1 1 1 1 1 1
469 %(output_displacement)i 0 0 0 0 0 0 %(output_rotation)i 0 0 0
471 '%(output_filename)s'
472 %(output_time_window)e
473 %(sfilter)s
474 %(output_slowness_min)e %(output_slowness_max)e
475 %(n_receiver_lines)i
476%(receiver_lines)s
477#-------------------------------------------------------------------------------------------
478#
479# LAYERED EARTH MODEL (IASP91)
480# ============================
481# 1. number of data lines of the layered model and selection for including
482# the physical dispersion according Kamamori & Anderson (1977)
483#-------------------------------------------------------------------------------------------
484 %(n_model_lines)i %(include_physical_dispersion)i
485#--------------------------------------------------------------------------------------------------------
486#
487# MODEL PARAMETERS
488# ================
489# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs
490#-------------------------------------------------------------------------------------------
491%(model_lines)s
492#---------------------------------end of all inputs-----------------------------------------
493''' # noqa
494 else:
496 template = '''# autogenerated QSSP input by qssp.py
497#
498# This is the input file of FORTRAN77 program "qssp2010" for calculating
499# synthetic seismograms of a self-gravitating, spherically symmetric,
500# isotropic and viscoelastic earth.
501#
502# by
503# Rongjiang Wang <wang@gfz-potsdam.de>
504# Helmholtz-Centre Potsdam
505# GFZ German Reseach Centre for Geosciences
506# Telegrafenberg, D-14473 Potsdam, Germany
507#
508# Last modified: Potsdam, July, 2010
509#
510# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
511# If not specified, SI Unit System is used overall!
512#
513# Coordinate systems:
514# spherical (r,t,p) with r = radial,
515# t = co-latitude,
516# p = east longitude.
517# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
518#
519# UNIFORM RECEIVER DEPTH
520# ======================
521# 1. uniform receiver depth [km]
522#-------------------------------------------------------------------------------------------
523 %(receiver_depth)e
524#-------------------------------------------------------------------------------------------
525#
526# TIME (FREQUENCY) SAMPLING
527# =========================
528# 1. time window [sec], sampling interval [sec]
529# 2. max. frequency [Hz] of Green's functions
530# 3. max. slowness [s/km] of Green's functions
531# Note: if the near-field static displacement is desired, the maximum slowness should not
532# be smaller than the S wave slowness in the receiver layer
533# 4. anti-aliasing factor (> 0 & < 1), if it is <= 0 or >= 1/e (~ 0.4), then
534# default value of 1/e is used (e.g., 0.1 = alias phases will be suppressed
535# to 10%% of their original amplitude)
536#
537# Note: The computation effort increases linearly the time window and
538# quadratically with the cut-off frequency.
539#-------------------------------------------------------------------------------------------
540 %(time_window)e %(sampling_interval)e
541 %(frequency_max)e
542 %(slowness_max)e
543 %(antialiasing_factor)e
544#-------------------------------------------------------------------------------------------
545#
546# SELF-GRAVITATING EFFECT
547# =======================
548# 1. the critical frequency [Hz] and the critical harmonic degree, below which
549# the self-gravitating effect should be included
550#-------------------------------------------------------------------------------------------
551 %(crit_frequency_sge)e %(crit_harmonic_degree_sge)i
552#-------------------------------------------------------------------------------------------
553#
554# WAVE TYPES
555# ==========
556# 1. selection (1/0 = yes/no) of speroidal modes (P-SV waves), selection of toroidal modes
557%(scutoff_doc)s
558#-------------------------------------------------------------------------------------------
559 %(spheroidal_modes)i %(toroidal_modes)i %(scutoff)s
560#-------------------------------------------------------------------------------------------
561# GREEN'S FUNCTION FILES
562# ======================
563# 1. number of discrete source depths, estimated radius of each source patch [km] and
564# directory for Green's functions
565# 2. list of the source depths [km], the respective file names of the Green's
566# functions (spectra) and the switch number (0/1) (0 = do not calculate
567# this Green's function because it exists already, 1 = calculate or update
568# this Green's function. Note: update is required if any of the above
569# parameters is changed)
570#-------------------------------------------------------------------------------------------
571 %(n_gf_lines)i %(source_patch_radius)e '%(gf_directory)s'
572 %(gf_lines)s
573#-------------------------------------------------------------------------------------------
574#
575# MULTI-EVENT SOURCE PARAMETERS
576# =============================
577# 1. number of discrete point sources and selection of the source data format
578# (1 or 2)
579# 2. list of the multi-event sources
580# Format 1:
581# M-Unit Mrr Mtt Mpp Mrt Mrp Mtp Lat Lon Depth T_origin T_rise
582# [Nm] [deg] [deg] [km] [sec] [sec]
583# Format 2:
584# Moment Strike Dip Rake Lat Lon Depth T_origin T_rise
585# [Nm] [deg] [deg] [deg] [deg] [deg] [km] [sec] [sec]
586#-------------------------------------------------------------------------------------------
587 %(n_source_lines)i %(point_source_type)i
588%(source_lines)s
589#-------------------------------------------------------------------------------------------
590#
591# RECEIVER PARAMETERS
592# ===================
593# 1. output file name and and selection of output format:
594# 1 = cartesian: vertical(z)/north(n)/east(e);
595# 2 = spherical: radial(r)/theta(t)/phi(p)
596# (Note: if output format 2 is selected, the epicenter (T_origin = 0)
597# 2. output time window [sec] (<= Green's function time window)
598%(sfilter_doc)s
599# 4. lower and upper slowness cut-off [s/km] (slowness band-pass filter)
600# 5. number of receiver
601# 6. list of the station parameters
602# Format:
603# Lat Lon Name Time_reduction
604# [deg] [deg] [sec]
605# (Note: Time_reduction = start time of the time window)
606#-------------------------------------------------------------------------------------------
607 '%(output_filename)s' %(output_format)i
608 %(output_time_window)e
609 %(sfilter)s
610 %(output_slowness_min)e %(output_slowness_max)e
611 %(n_receiver_lines)i
612%(receiver_lines)s
613#-------------------------------------------------------------------------------------------
614#
615# LAYERED EARTH MODEL (IASP91)
616# ============================
617# 1. number of data lines of the layered model and selection for including
618# the physical dispersion according Kamamori & Anderson (1977)
619#-------------------------------------------------------------------------------------------
620 %(n_model_lines)i %(include_physical_dispersion)i
621#-------------------------------------------------------------------------------------------
622#
623# MULTILAYERED MODEL PARAMETERS (source site)
624# ===========================================
625# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs
626#-------------------------------------------------------------------------------------------
627%(model_lines)s
628#---------------------------------end of all inputs-----------------------------------------
629''' # noqa
631 return (template % d).encode('ascii')
634class QSSPError(gf.store.StoreError):
635 pass
638class Interrupted(gf.store.StoreError):
639 def __str__(self):
640 return 'Interrupted.'
643class QSSPRunner(object):
645 def __init__(self, tmp=None, keep_tmp=False):
647 self.tempdir = mkdtemp(prefix='qssprun-', dir=tmp)
648 self.keep_tmp = keep_tmp
649 self.config = None
651 def run(self, config):
652 self.config = config
654 input_fn = pjoin(self.tempdir, 'input')
656 with open(input_fn, 'wb') as f:
657 input_str = config.string_for_config()
658 logger.debug('===== begin qssp input =====\n'
659 '%s===== end qssp input =====' % input_str.decode())
660 f.write(input_str)
662 program = program_bins['qssp.%s' % config.qssp_version]
664 old_wd = os.getcwd()
666 os.chdir(self.tempdir)
668 interrupted = []
670 def signal_handler(signum, frame):
671 os.kill(proc.pid, signal.SIGTERM)
672 interrupted.append(True)
674 original = signal.signal(signal.SIGINT, signal_handler)
675 try:
676 try:
677 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
678 except OSError:
679 os.chdir(old_wd)
680 raise QSSPError(
681 '''could not start qssp executable: "%s"
682Available fomosto backends and download links to the modelling codes are listed
683on
685 https://pyrocko.org/docs/current/apps/fomosto/backends.html
687''' % program)
689 (output_str, error_str) = proc.communicate(b'input\n')
691 finally:
692 signal.signal(signal.SIGINT, original)
694 if interrupted:
695 raise KeyboardInterrupt()
697 logger.debug('===== begin qssp output =====\n'
698 '%s===== end qssp output =====' % output_str.decode())
700 errmess = []
701 if proc.returncode != 0:
702 errmess.append(
703 'qssp had a non-zero exit state: %i' % proc.returncode)
704 if error_str:
706 logger.warning(
707 'qssp emitted something via stderr: \n\n%s'
708 % error_str.decode())
710 # errmess.append('qssp emitted something via stderr')
711 if output_str.lower().find(b'error') != -1:
712 errmess.append("the string 'error' appeared in qssp output")
714 if errmess:
715 os.chdir(old_wd)
716 raise QSSPError('''
717===== begin qssp input =====
718%s===== end qssp input =====
719===== begin qssp output =====
720%s===== end qssp output =====
721===== begin qssp error =====
722%s===== end qssp error =====
723%s
724qssp has been invoked as "%s"'''.lstrip() % (
725 input_str.decode(),
726 output_str.decode(),
727 error_str.decode(),
728 '\n'.join(errmess),
729 program))
731 self.qssp_output = output_str
732 self.qssp_error = error_str
734 os.chdir(old_wd)
736 def get_traces(self):
738 fns = self.config.get_output_filenames(self.tempdir)
739 traces = {}
740 for comp, fn in zip(self.config.components, fns):
741 data = num.loadtxt(fn, skiprows=1, dtype=float)
742 nsamples, ntraces = data.shape
743 ntraces -= 1
744 deltat = (data[-1, 0] - data[0, 0])/(nsamples-1)
745 toffset = data[0, 0]
746 for itrace in range(ntraces):
747 rec = self.config.receivers[itrace]
748 tmin = rec.tstart + toffset
749 tr = trace.Trace(
750 '', '%04i' % itrace, '', comp,
751 tmin=tmin, deltat=deltat, ydata=data[:, itrace+1],
752 meta=dict(distance=rec.distance))
754 traces[itrace, comp] = tr
756 if self.config.qssp_version == 'ppeg2017':
757 for itrace in range(ntraces):
758 for c in 'nez':
759 tr_accel = traces[itrace, '_acce_' + c]
760 tr_gravi = traces[itrace, '_gravitation_' + c]
761 tr_ag = tr_accel.copy()
762 tr_ag.ydata -= tr_gravi.ydata
763 tr_ag.set_codes(channel='ag_' + c)
764 tr_ag.meta = tr_accel.meta
765 traces[itrace, 'ag_' + c] = tr_ag
767 traces = list(traces.values())
768 traces.sort(key=lambda tr: tr.nslc_id)
769 return traces
771 def __del__(self):
772 if self.tempdir:
773 if not self.keep_tmp:
774 shutil.rmtree(self.tempdir)
775 self.tempdir = None
776 else:
777 logger.warning(
778 'not removing temporary directory: %s' % self.tempdir)
781class QSSPGFBuilder(gf.builder.Builder):
782 nsteps = 2
784 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
785 force=False):
787 self.store = gf.store.Store(store_dir, 'w')
788 baseconf = self.store.get_extra('qssp')
789 if baseconf.qssp_version in ('2017', '2020'):
790 if self.store.config.stored_quantity == "rotation":
791 self.gfmapping = [
792 (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
793 {'_rota_n': (0, -1), '_rota_e': (3, -1),
794 '_rota_z': (5, -1)}),
795 (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
796 {'_rota_n': (1, -1), '_rota_e': (4, -1),
797 '_rota_z': (6, -1)}),
798 (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
799 {'_rota_n': (2, -1), '_rota_z': (7, -1)}),
800 (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
801 {'_rota_n': (8, -1), '_rota_z': (9, -1)}),
802 ]
803 else:
804 self.gfmapping = [
805 (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
806 {'_disp_n': (0, -1), '_disp_e': (3, -1),
807 '_disp_z': (5, -1)}),
808 (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
809 {'_disp_n': (1, -1), '_disp_e': (4, -1),
810 '_disp_z': (6, -1)}),
811 (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
812 {'_disp_n': (2, -1), '_disp_z': (7, -1)}),
813 (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
814 {'_disp_n': (8, -1), '_disp_z': (9, -1)}),
815 ]
817 elif baseconf.qssp_version == 'ppeg2017':
818 self.gfmapping = [
819 (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
820 {'ag_n': (0, -1), 'ag_e': (3, -1), 'ag_z': (5, -1)}),
821 (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
822 {'ag_n': (1, -1), 'ag_e': (4, -1), 'ag_z': (6, -1)}),
823 (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
824 {'ag_n': (2, -1), 'ag_z': (7, -1)}),
825 (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
826 {'ag_n': (8, -1), 'ag_z': (9, -1)}),
827 ]
828 else:
829 self.gfmapping = [
830 (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
831 {'un': (0, -1), 'ue': (3, -1), 'uz': (5, -1)}),
832 (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
833 {'un': (1, -1), 'ue': (4, -1), 'uz': (6, -1)}),
834 (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
835 {'un': (2, -1), 'uz': (7, -1)}),
836 (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
837 {'un': (8, -1), 'uz': (9, -1)}),
838 ]
840 if step == 0:
841 block_size = (1, 1, self.store.config.ndistances)
842 else:
843 if block_size is None:
844 block_size = (1, 1, 51)
846 if len(self.store.config.ns) == 2:
847 block_size = block_size[1:]
849 gf.builder.Builder.__init__(
850 self, self.store.config, step, block_size=block_size, force=force)
852 conf = QSSPConfigFull(**baseconf.items())
853 conf.gf_directory = pjoin(store_dir, 'qssp_green')
854 conf.earthmodel_1d = self.store.config.earthmodel_1d
855 deltat = self.store.config.deltat
856 if 'time_window' not in shared:
857 d = self.store.make_timing_params(
858 conf.time_region[0], conf.time_region[1],
859 force=force)
861 tmax = math.ceil(d['tmax'] / deltat) * deltat
862 tmin = math.floor(d['tmin'] / deltat) * deltat
864 shared['time_window'] = tmax - tmin
865 shared['tstart'] = tmin
867 self.tstart = shared['tstart']
868 conf.time_window = shared['time_window']
870 self.tmp = tmp
871 if self.tmp is not None:
872 util.ensuredir(self.tmp)
874 util.ensuredir(conf.gf_directory)
876 self.qssp_config = conf
878 def cleanup(self):
879 self.store.close()
881 def work_block(self, iblock):
882 if len(self.store.config.ns) == 2:
883 (sz, firstx), (sz, lastx), (ns, nx) = \
884 self.get_block_extents(iblock)
886 rz = self.store.config.receiver_depth
887 else:
888 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
889 self.get_block_extents(iblock)
891 gf_filename = 'GF_%gkm_%gkm' % (sz/km, rz/km)
893 conf = copy.deepcopy(self.qssp_config)
895 gf_path = os.path.join(conf.gf_directory, '?_' + gf_filename)
897 if self.step == 0 and len(glob.glob(gf_path)) > 0:
898 logger.info(
899 'Skipping step %i / %i, block %i / %i (GF already exists)'
900 % (self.step+1, self.nsteps, iblock+1, self.nblocks))
902 return
904 logger.info(
905 'Starting step %i / %i, block %i / %i' %
906 (self.step+1, self.nsteps, iblock+1, self.nblocks))
908 tbeg = time.time()
910 runner = QSSPRunner(tmp=self.tmp)
912 conf.receiver_depth = rz/km
913 conf.stored_quantity = self.store.config.stored_quantity
914 conf.sampling_interval = 1.0 / self.gf_config.sample_rate
915 dx = self.gf_config.distance_delta
917 if self.step == 0:
918 distances = [firstx]
919 else:
920 distances = num.linspace(firstx, firstx + (nx-1)*dx, nx)
922 conf.receivers = [
923 QSSPReceiver(
924 lat=90-d*cake.m2d,
925 lon=180.,
926 tstart=self.tstart,
927 distance=d)
929 for d in distances]
931 if self.step == 0:
932 gf_filename = 'TEMP' + gf_filename[2:]
934 gfs = [QSSPGreen(
935 filename=gf_filename,
936 depth=sz/km,
937 calculate=(self.step == 0))]
939 conf.greens_functions = gfs
941 trise = 0.001*conf.sampling_interval # make it short (delta impulse)
943 if self.step == 0:
944 conf.sources = [QSSPSourceMT(
945 lat=90-0.001*dx*cake.m2d,
946 lon=0.0,
947 trise=trise,
948 torigin=0.0)]
950 runner.run(conf)
951 gf_path = os.path.join(conf.gf_directory, '?_' + gf_filename)
952 for s in glob.glob(gf_path):
953 d = s.replace('TEMP_', 'GF_')
954 os.rename(s, d)
956 else:
957 for mt, gfmap in self.gfmapping[
958 :[3, 4][self.gf_config.ncomponents == 10]]:
959 m = mt.m_up_south_east()
961 conf.sources = [QSSPSourceMT(
962 lat=90-0.001*dx*cake.m2d,
963 lon=0.0,
964 mrr=m[0, 0], mtt=m[1, 1], mpp=m[2, 2],
965 mrt=m[0, 1], mrp=m[0, 2], mtp=m[1, 2],
966 trise=trise,
967 torigin=0.0)]
969 runner.run(conf)
971 rawtraces = runner.get_traces()
973 interrupted = []
975 def signal_handler(signum, frame):
976 interrupted.append(True)
978 original = signal.signal(signal.SIGINT, signal_handler)
979 self.store.lock()
980 duplicate_inserts = 0
981 try:
982 for itr, tr in enumerate(rawtraces):
983 if tr.channel in gfmap:
985 x = tr.meta['distance']
986 ig, factor = gfmap[tr.channel]
988 if len(self.store.config.ns) == 2:
989 args = (sz, x, ig)
990 else:
991 args = (rz, sz, x, ig)
993 if conf.cut:
994 tmin = self.store.t(conf.cut[0], args[:-1])
995 tmax = self.store.t(conf.cut[1], args[:-1])
996 if None in (tmin, tmax):
997 continue
999 tr.chop(tmin, tmax)
1001 tmin = tr.tmin
1002 tmax = tr.tmax
1004 if conf.fade:
1005 ta, tb, tc, td = [
1006 self.store.t(v, args[:-1])
1007 for v in conf.fade]
1009 if None in (ta, tb, tc, td):
1010 continue
1012 if not (ta <= tb and tb <= tc and tc <= td):
1013 raise QSSPError(
1014 'invalid fade configuration')
1016 t = tr.get_xdata()
1017 fin = num.interp(t, [ta, tb], [0., 1.])
1018 fout = num.interp(t, [tc, td], [1., 0.])
1019 anti_fin = 1. - fin
1020 anti_fout = 1. - fout
1022 y = tr.ydata
1024 sum_anti_fin = num.sum(anti_fin)
1025 sum_anti_fout = num.sum(anti_fout)
1027 if conf.nonzero_fade_in \
1028 and sum_anti_fin != 0.0:
1029 yin = num.sum(anti_fin*y) / sum_anti_fin
1030 else:
1031 yin = 0.0
1033 if conf.nonzero_fade_out \
1034 and sum_anti_fout != 0.0:
1035 yout = num.sum(anti_fout*y) / sum_anti_fout
1036 else:
1037 yout = 0.0
1039 y2 = anti_fin*yin + fin*fout*y + anti_fout*yout
1041 if conf.relevel_with_fade_in:
1042 y2 -= yin
1044 tr.set_ydata(y2)
1046 gf_tr = gf.store.GFTrace.from_trace(tr)
1047 gf_tr.data *= factor
1049 try:
1050 self.store.put(args, gf_tr)
1051 except gf.store.DuplicateInsert:
1052 duplicate_inserts += 1
1054 finally:
1055 if duplicate_inserts:
1056 logger.warning(
1057 '%i insertions skipped (duplicates)'
1058 % duplicate_inserts)
1060 self.store.unlock()
1061 signal.signal(signal.SIGINT, original)
1063 if interrupted:
1064 raise KeyboardInterrupt()
1066 tend = time.time()
1067 logger.info(
1068 'Done with step %i / %i, block %i / %i, wallclock time: %.0f s' % (
1069 self.step+1, self.nsteps, iblock+1, self.nblocks, tend-tbeg))
1072km = 1000.
1075def init(store_dir, variant, config_params=None):
1076 if variant is None:
1077 variant = '2010beta'
1079 if ('qssp.' + variant) not in program_bins:
1080 raise gf.store.StoreError('unsupported qssp variant: %s' % variant)
1082 qssp = QSSPConfig(qssp_version=variant)
1083 if variant != 'ppeg2017':
1084 qssp.time_region = (
1085 gf.Timing('begin-50'),
1086 gf.Timing('end+100'))
1088 qssp.cut = (
1089 gf.Timing('begin-50'),
1090 gf.Timing('end+100'))
1092 else: # variant == 'ppeg2017':
1093 qssp.frequency_max = 0.5
1094 qssp.time_region = [
1095 gf.Timing('-100'), gf.Timing('{stored:begin}+100')]
1096 qssp.cut = [
1097 gf.Timing('-100'), gf.Timing('{stored:begin}+100')]
1098 qssp.antialiasing_factor = 1.0e-10
1099 qssp.toroidal_modes = False
1100 qssp.cutoff_harmonic_degree_min = 2500
1101 qssp.cutoff_harmonic_degree_max = 2500
1102 qssp.crit_frequency_sge = 5.0
1103 qssp.crit_harmonic_degree_sge = 50000
1104 qssp.source_patch_radius = 10.0
1105 qssp.bandpass_order = 6
1106 qssp.bandpass_corner_low = 0.0
1107 qssp.bandpass_corner_high = 0.125
1109 store_id = os.path.basename(os.path.realpath(store_dir))
1110 if variant == 'ppeg2017':
1111 quantity = 'acceleration'
1112 else:
1113 quantity = None
1115 if variant == 'ppeg2017':
1116 sample_rate = 4.0
1117 else:
1118 sample_rate = 0.2
1120 d = dict(
1121 id=store_id,
1122 ncomponents=10,
1123 component_scheme='elastic10',
1124 stored_quantity=quantity,
1125 sample_rate=sample_rate,
1126 receiver_depth=0*km,
1127 source_depth_min=10*km,
1128 source_depth_max=20*km,
1129 source_depth_delta=10*km,
1130 distance_min=100*km,
1131 distance_max=1000*km,
1132 distance_delta=10*km,
1133 earthmodel_1d=cake.load_model(),
1134 modelling_code_id='qssp',
1135 tabulated_phases=[
1136 gf.meta.TPDef(
1137 id='begin',
1138 definition='p,P,p\\,P\\,Pv_(cmb)p'),
1139 gf.meta.TPDef(
1140 id='end',
1141 definition='2.5'),
1142 gf.meta.TPDef(
1143 id='P',
1144 definition='!P'),
1145 gf.meta.TPDef(
1146 id='S',
1147 definition='!S'),
1148 gf.meta.TPDef(
1149 id='p',
1150 definition='!p'),
1151 gf.meta.TPDef(
1152 id='s',
1153 definition='!s')])
1155 if config_params is not None:
1156 d.update(config_params)
1158 config = gf.meta.ConfigTypeA(**d)
1160 config.validate()
1161 return gf.store.Store.create_editables(
1162 store_dir,
1163 config=config,
1164 extra={'qssp': qssp})
1167def build(
1168 store_dir,
1169 force=False,
1170 nworkers=None,
1171 continue_=False,
1172 step=None,
1173 iblock=None):
1175 return QSSPGFBuilder.build(
1176 store_dir, force=force, nworkers=nworkers, continue_=continue_,
1177 step=step, iblock=iblock)
1180def get_conf(
1181 store_dir,
1182 force=False,
1183 nworkers=None,
1184 continue_=False,
1185 step=None,
1186 iblock=None):
1188 return QSSPGFBuilder.get_conf()