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 math
11import signal
13from pyrocko import trace, cake, gf
14from pyrocko.ahfullgreen import add_seismogram, Impulse
15from pyrocko.moment_tensor import MomentTensor, symmat6
17km = 1000.
19guts_prefix = 'pf'
21logger = logging.getLogger('pyrocko.fomosto.ahfullgreen')
23# how to call the programs
24program_bins = {
25 'ahfullgreen': 'ahfullgreen',
26}
28components = 'r t z'.split()
31def example_model():
32 material = cake.Material(vp=3000., vs=1000., rho=3000., qp=200., qs=100.)
33 layer = cake.HomogeneousLayer(
34 ztop=0., zbot=30*km, m=material, name='fullspace')
35 mod = cake.LayeredModel()
36 mod.append(layer)
37 return mod
40class AhfullgreenError(gf.store.StoreError):
41 pass
44def make_traces(material, source_mech, deltat, norths, easts,
45 source_depth, receiver_depth):
47 if isinstance(source_mech, MomentTensor):
48 m6 = source_mech.m6()
49 f = (0., 0., 0.)
50 elif isinstance(source_mech, tuple):
51 m6 = (0., 0., 0., 0., 0., 0.)
52 f = source_mech
54 npad = 120
56 traces = []
57 for i_distance, (north, east) in enumerate(zip(norths, easts)):
58 d3d = math.sqrt(
59 north**2 + east**2 + (receiver_depth - source_depth)**2)
61 tmin = (math.floor(d3d / material.vp / deltat) - npad) * deltat
62 tmax = (math.ceil(d3d / material.vs / deltat) + npad) * deltat
63 ns = int(round((tmax - tmin) / deltat))
65 outx = num.zeros(ns)
66 outy = num.zeros(ns)
67 outz = num.zeros(ns)
69 x = (north, east, receiver_depth-source_depth)
71 add_seismogram(
72 material.vp, material.vs, material.rho, material.qp, material.qs,
73 x, f, m6, 'displacement',
74 deltat, tmin, outx, outy, outz,
75 stf=Impulse())
77 for i_comp, o in enumerate((outx, outy, outz)):
78 comp = components[i_comp]
79 tr = trace.Trace('', '%04i' % i_distance, '', comp,
80 tmin=tmin, ydata=o, deltat=deltat,
81 meta=dict(
82 distance=math.sqrt(north**2 + east**2),
83 azimuth=0.))
85 traces.append(tr)
87 return traces
90class AhfullGFBuilder(gf.builder.Builder):
91 def __init__(self, store_dir, step, shared, block_size=None, force=False):
93 self.store = gf.store.Store(store_dir, 'w')
95 if block_size is None:
96 block_size = (1, 1, 2000)
98 if len(self.store.config.ns) == 2:
99 block_size = block_size[1:]
101 gf.builder.Builder.__init__(
102 self, self.store.config, step, block_size=block_size, force=force)
104 def cleanup(self):
105 self.store.close()
107 def work_block(self, index):
108 if len(self.store.config.ns) == 2:
109 (sz, firstx), (sz, lastx), (ns, nx) = \
110 self.get_block_extents(index)
112 rz = self.store.config.receiver_depth
113 else:
114 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
115 self.get_block_extents(index)
117 logger.info('Starting block %i / %i' %
118 (index+1, self.nblocks))
120 dx = self.gf_config.distance_delta
122 distances = num.linspace(firstx, firstx + (nx-1)*dx, nx).tolist()
124 mmt1 = (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
125 {'r': (0, +1), 't': (3, +1), 'z': (5, +1)})
126 mmt2 = (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
127 {'r': (1, +1), 't': (4, +1), 'z': (6, +1)})
128 mmt3 = (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
129 {'r': (2, +1), 'z': (7, +1)})
130 mmt4 = (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
131 {'r': (8, +1), 'z': (9, +1)})
132 mmt0 = (MomentTensor(m=symmat6(1, 1, 1, 0, 0, 0)),
133 {'r': (0, +1), 'z': (1, +1)})
135 component_scheme = self.store.config.component_scheme
137 if component_scheme == 'elastic8':
138 gfmapping = [mmt1, mmt2, mmt3]
140 elif component_scheme == 'elastic10':
141 gfmapping = [mmt1, mmt2, mmt3, mmt4]
143 elif component_scheme == 'elastic2':
144 gfmapping = [mmt0]
146 elif component_scheme == 'elastic5':
147 gfmapping = [
148 ((1., 1., 0.), {'r': (1, +1), 'z': (4, +1), 't': (2, +1)}),
149 ((0., 0., 1.), {'r': (0, +1), 'z': (3, +1)})]
151 else:
152 raise gf.UnavailableScheme(
153 'fomosto backend "ahfullgreen" cannot handle component scheme '
154 '"%s"' % component_scheme)
156 for source_mech, gfmap in gfmapping:
158 rawtraces = make_traces(
159 self.store.config.earthmodel_1d.require_homogeneous(),
160 source_mech, 1.0/self.store.config.sample_rate,
161 distances, num.zeros_like(distances), sz, rz)
163 interrupted = []
165 def signal_handler(signum, frame):
166 interrupted.append(True)
168 original = signal.signal(signal.SIGINT, signal_handler)
169 self.store.lock()
170 duplicate_inserts = 0
171 try:
172 for itr, tr in enumerate(rawtraces):
173 if tr.channel not in gfmap:
174 logger.debug('%s not in gfmap' % tr.channel)
175 continue
177 x = tr.meta['distance']
178 if x > firstx + (nx-1)*dx:
179 logger.error("x out of range")
180 continue
182 ig, factor = gfmap[tr.channel]
184 if len(self.store.config.ns) == 2:
185 args = (sz, x, ig)
186 else:
187 args = (rz, sz, x, ig)
189 tr = tr.snap()
191 gf_tr = gf.store.GFTrace.from_trace(tr)
192 gf_tr.data *= factor
194 try:
195 self.store.put(args, gf_tr)
196 except gf.store.DuplicateInsert:
197 duplicate_inserts += 1
199 finally:
200 if duplicate_inserts:
201 logger.warn('%i insertions skipped (duplicates)' %
202 duplicate_inserts)
204 self.store.unlock()
205 signal.signal(signal.SIGINT, original)
207 if interrupted:
208 raise KeyboardInterrupt()
210 logger.info('Done with block %i / %i' %
211 (index+1, self.nblocks))
214def init(store_dir, variant):
215 assert variant is None
217 modelling_code_id = 'ahfullgreen'
219 store_id = os.path.basename(os.path.realpath(store_dir))
221 config = gf.meta.ConfigTypeA(
222 id=store_id,
223 ncomponents=10,
224 sample_rate=20.,
225 receiver_depth=0*km,
226 source_depth_min=1*km,
227 source_depth_max=10*km,
228 source_depth_delta=1*km,
229 distance_min=1*km,
230 distance_max=20*km,
231 distance_delta=1*km,
232 earthmodel_1d=example_model(),
233 modelling_code_id=modelling_code_id,
234 tabulated_phases=[
235 gf.meta.TPDef(
236 id='anyP',
237 definition='P,p,\\P,\\p'),
238 gf.meta.TPDef(
239 id='anyS',
240 definition='S,s,\\S,\\s')])
242 config.validate()
243 return gf.store.Store.create_editables(
244 store_dir, config=config)
247def build(store_dir, force=False, nworkers=None, continue_=False, step=None,
248 iblock=None):
250 return AhfullGFBuilder.build(
251 store_dir, force=force, nworkers=nworkers, continue_=continue_,
252 step=step, iblock=iblock)