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