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.guts import Object, clone
14from pyrocko import trace, cake, gf
15from pyrocko.ahfullgreen import add_seismogram, AhfullgreenSTFImpulse, \
16 AhfullgreenSTF
17from pyrocko.moment_tensor import MomentTensor, symmat6
19km = 1000.
21guts_prefix = 'pf'
23logger = logging.getLogger('pyrocko.fomosto.ahfullgreen')
25# how to call the programs
26program_bins = {
27 'ahfullgreen': 'ahfullgreen',
28}
30components = 'x y z'.split()
33def example_model():
34 material = cake.Material(vp=3000., vs=1000., rho=3000., qp=200., qs=100.)
35 layer = cake.HomogeneousLayer(
36 ztop=0., zbot=30*km, m=material, name='fullspace')
37 mod = cake.LayeredModel()
38 mod.append(layer)
39 return mod
42class AhfullgreenError(gf.store.StoreError):
43 pass
46class AhfullgreenConfig(Object):
47 stf = AhfullgreenSTF.T(default=AhfullgreenSTFImpulse.D())
50def make_traces(material, source_mech, deltat, norths, easts,
51 source_depth, receiver_depth, stf):
53 if isinstance(source_mech, MomentTensor):
54 m6 = source_mech.m6()
55 f = (0., 0., 0.)
56 elif isinstance(source_mech, tuple):
57 m6 = (0., 0., 0., 0., 0., 0.)
58 f = source_mech
60 npad = 120
62 traces = []
63 for isource, (north, east) in enumerate(zip(norths, easts)):
64 d3d = math.sqrt(
65 north**2 + east**2 + (receiver_depth - source_depth)**2)
67 tmin = (math.floor(d3d / material.vp / deltat) - npad) * deltat
68 tmax = (math.ceil(d3d / material.vs / deltat) + npad) * deltat
69 ns = int(round((tmax - tmin) / deltat))
71 outx = num.zeros(ns)
72 outy = num.zeros(ns)
73 outz = num.zeros(ns)
75 x = (north, east, receiver_depth-source_depth)
77 add_seismogram(
78 material.vp, material.vs, material.rho, material.qp, material.qs,
79 x, f, m6, 'displacement',
80 deltat, tmin, outx, outy, outz,
81 stf=stf)
83 for i_comp, o in enumerate((outx, outy, outz)):
84 comp = components[i_comp]
85 tr = trace.Trace('', '%04i' % isource, '', comp,
86 tmin=tmin, ydata=o, deltat=deltat,
87 meta=dict(isource=isource))
89 traces.append(tr)
91 return traces
94class AhfullGFBuilder(gf.builder.Builder):
95 def __init__(self, store_dir, step, shared, block_size=None, force=False):
97 self.store = gf.store.Store(store_dir, 'w')
99 block_size = {
100 'A': (1, 2000),
101 'B': (1, 1, 2000),
102 'C': (1, 1, 2)}[self.store.config.short_type]
104 gf.builder.Builder.__init__(
105 self, self.store.config, step, block_size=block_size, force=force)
107 self.ahfullgreen_config = self.store.get_extra('ahfullgreen')
109 def cleanup(self):
110 self.store.close()
112 def work_block(self, index):
113 store_type = self.store.config.short_type
115 if store_type == 'A':
116 (sz, firstx), (sz, lastx), (ns, nx) = \
117 self.get_block_extents(index)
119 rz = self.store.config.receiver_depth
120 sy = 0.0
121 elif store_type == 'B':
122 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
123 self.get_block_extents(index)
124 sy = 0.0
125 elif store_type == 'C':
126 (sz, sy, firstx), (sz, sy, lastx), (nz, ny, nx) = \
127 self.get_block_extents(index)
128 rz = self.store.config.receiver.depth
130 logger.info('Starting block %i / %i' %
131 (index+1, self.nblocks))
133 dx = self.gf_config.deltas[-1]
135 xs = num.linspace(firstx, firstx + (nx-1)*dx, nx).tolist()
137 mmt1 = (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
138 {'x': (0, +1), 'y': (3, +1), 'z': (5, +1)})
139 mmt2 = (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
140 {'x': (1, +1), 'y': (4, +1), 'z': (6, +1)})
141 mmt3 = (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
142 {'x': (2, +1), 'z': (7, +1)})
143 mmt4 = (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
144 {'x': (8, +1), 'z': (9, +1)})
145 mmt0 = (MomentTensor(m=symmat6(1, 1, 1, 0, 0, 0)),
146 {'x': (0, +1), 'z': (1, +1)})
148 component_scheme = self.store.config.component_scheme
150 if component_scheme == 'elastic8':
151 gfmapping = [mmt1, mmt2, mmt3]
153 elif component_scheme == 'elastic10':
154 gfmapping = [mmt1, mmt2, mmt3, mmt4]
156 elif component_scheme == 'elastic2':
157 gfmapping = [mmt0]
159 elif component_scheme == 'elastic5':
160 gfmapping = [
161 ((1., 1., 0.), {'x': (1, +1), 'z': (4, +1), 'y': (2, +1)}),
162 ((0., 0., 1.), {'x': (0, +1), 'z': (3, +1)})]
163 elif component_scheme == 'elastic18':
164 gfmapping = []
165 for im in range(6):
166 m6 = [0.0] * 6
167 m6[im] = 1.0
169 gfmapping.append(
170 (MomentTensor(m=symmat6(*m6)),
171 {'x': (im, +1),
172 'y': (6+im, +1),
173 'z': (12+im, +1)}))
175 else:
176 raise gf.UnavailableScheme(
177 'fomosto backend "ahfullgreen" cannot handle component scheme '
178 '"%s"' % component_scheme)
180 for source_mech, gfmap in gfmapping:
182 if component_scheme != 'elastic18':
183 norths = xs
184 easts = num.zeros_like(xs)
185 else:
186 receiver = self.store.config.receiver
187 data = []
188 for x in xs:
189 source = clone(self.store.config.source_origin)
190 source.north_shift = x
191 source.east_shift = sy
192 source.depth = sz
193 data.append(receiver.offset_to(source))
195 norths, easts = -num.array(data).T
197 rawtraces = make_traces(
198 self.store.config.earthmodel_1d.require_homogeneous(),
199 source_mech, 1.0/self.store.config.sample_rate,
200 norths, easts, sz, rz, self.ahfullgreen_config.stf)
202 interrupted = []
204 def signal_handler(signum, frame):
205 interrupted.append(True)
207 original = signal.signal(signal.SIGINT, signal_handler)
208 self.store.lock()
209 duplicate_inserts = 0
210 try:
211 for itr, tr in enumerate(rawtraces):
212 if tr.channel not in gfmap:
213 logger.debug('%s not in gfmap' % tr.channel)
214 continue
216 x = xs[tr.meta['isource']]
217 if x > firstx + (nx-1)*dx:
218 logger.error("x out of range")
219 continue
221 ig, factor = gfmap[tr.channel]
223 args = {
224 'A': (sz, x, ig),
225 'B': (rz, sz, x, ig),
226 'C': (sz, sy, x, ig)}[store_type]
228 tr = tr.snap()
230 gf_tr = gf.store.GFTrace.from_trace(tr)
231 gf_tr.data *= factor
233 try:
234 self.store.put(args, gf_tr)
235 except gf.store.DuplicateInsert:
236 duplicate_inserts += 1
238 finally:
239 if duplicate_inserts:
240 logger.warning('%i insertions skipped (duplicates)' %
241 duplicate_inserts)
243 self.store.unlock()
244 signal.signal(signal.SIGINT, original)
246 if interrupted:
247 raise KeyboardInterrupt()
249 logger.info('Done with block %i / %i' %
250 (index+1, self.nblocks))
253def init(store_dir, variant, config_params=None):
254 assert variant is None
256 modelling_code_id = 'ahfullgreen'
258 store_id = os.path.basename(os.path.realpath(store_dir))
260 ahfullgreen = AhfullgreenConfig()
262 d = dict(
263 id=store_id,
264 ncomponents=10,
265 sample_rate=20.,
266 receiver_depth=0*km,
267 source_depth_min=1*km,
268 source_depth_max=10*km,
269 source_depth_delta=1*km,
270 distance_min=1*km,
271 distance_max=20*km,
272 distance_delta=1*km,
273 earthmodel_1d=example_model(),
274 modelling_code_id=modelling_code_id,
275 tabulated_phases=[
276 gf.meta.TPDef(
277 id='anyP',
278 definition='P,p,\\P,\\p'),
279 gf.meta.TPDef(
280 id='anyS',
281 definition='S,s,\\S,\\s')])
283 if config_params is not None:
284 d.update(config_params)
286 config = gf.meta.ConfigTypeA(**d)
287 config.validate()
288 return gf.store.Store.create_editables(
289 store_dir, config=config, extra={'ahfullgreen': ahfullgreen})
292def build(store_dir, force=False, nworkers=None, continue_=False, step=None,
293 iblock=None):
295 return AhfullGFBuilder.build(
296 store_dir, force=force, nworkers=nworkers, continue_=continue_,
297 step=step, iblock=iblock)