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 math
10import signal
12from pyrocko.guts import Object, clone
13from pyrocko import trace, cake, gf
14from pyrocko.ahfullgreen import add_seismogram, AhfullgreenSTFImpulse, \
15 AhfullgreenSTF
16from pyrocko.moment_tensor import MomentTensor, symmat6
18km = 1000.
20guts_prefix = 'pf'
22logger = logging.getLogger('pyrocko.fomosto.ahfullgreen')
24# how to call the programs
25program_bins = {
26 'ahfullgreen': 'ahfullgreen',
27}
29components = 'x y z'.split()
32def example_model():
33 material = cake.Material(vp=3000., vs=1000., rho=3000., qp=200., qs=100.)
34 layer = cake.HomogeneousLayer(
35 ztop=0., zbot=30*km, m=material, name='fullspace')
36 mod = cake.LayeredModel()
37 mod.append(layer)
38 return mod
41class AhfullgreenError(gf.store.StoreError):
42 pass
45class AhfullgreenConfig(Object):
46 stf = AhfullgreenSTF.T(default=AhfullgreenSTFImpulse.D())
49def make_traces(material, source_mech, deltat, norths, easts,
50 source_depth, receiver_depth, stf):
52 if isinstance(source_mech, MomentTensor):
53 m6 = source_mech.m6()
54 f = (0., 0., 0.)
55 elif isinstance(source_mech, tuple):
56 m6 = (0., 0., 0., 0., 0., 0.)
57 f = source_mech
59 npad = 120
61 traces = []
62 for isource, (north, east) in enumerate(zip(norths, easts)):
63 d3d = math.sqrt(
64 north**2 + east**2 + (receiver_depth - source_depth)**2)
66 tmin = (math.floor(d3d / material.vp / deltat) - npad) * deltat
67 tmax = (math.ceil(d3d / material.vs / deltat) + npad) * deltat
68 ns = int(round((tmax - tmin) / deltat))
70 outx = num.zeros(ns)
71 outy = num.zeros(ns)
72 outz = num.zeros(ns)
74 x = (north, east, receiver_depth-source_depth)
76 add_seismogram(
77 material.vp, material.vs, material.rho, material.qp, material.qs,
78 x, f, m6, 'displacement',
79 deltat, tmin, outx, outy, outz,
80 stf=stf)
82 for i_comp, o in enumerate((outx, outy, outz)):
83 comp = components[i_comp]
84 tr = trace.Trace('', '%04i' % isource, '', comp,
85 tmin=tmin, ydata=o, deltat=deltat,
86 meta=dict(isource=isource))
88 traces.append(tr)
90 return traces
93class AhfullGFBuilder(gf.builder.Builder):
94 def __init__(self, store_dir, step, shared, block_size=None, force=False):
96 self.store = gf.store.Store(store_dir, 'w')
98 block_size = {
99 'A': (1, 2000),
100 'B': (1, 1, 2000),
101 'C': (1, 1, 2)}[self.store.config.short_type]
103 gf.builder.Builder.__init__(
104 self, self.store.config, step, block_size=block_size, force=force)
106 self.ahfullgreen_config = self.store.get_extra('ahfullgreen')
108 def cleanup(self):
109 self.store.close()
111 def work_block(self, index):
112 store_type = self.store.config.short_type
114 if store_type == 'A':
115 (sz, firstx), (sz, lastx), (ns, nx) = \
116 self.get_block_extents(index)
118 rz = self.store.config.receiver_depth
119 sy = 0.0
120 elif store_type == 'B':
121 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
122 self.get_block_extents(index)
123 sy = 0.0
124 elif store_type == 'C':
125 (sz, sy, firstx), (sz, sy, lastx), (nz, ny, nx) = \
126 self.get_block_extents(index)
127 rz = self.store.config.receiver.depth
129 logger.info('Starting block %i / %i' %
130 (index+1, self.nblocks))
132 dx = self.gf_config.deltas[-1]
134 xs = num.linspace(firstx, firstx + (nx-1)*dx, nx).tolist()
136 mmt1 = (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
137 {'x': (0, +1), 'y': (3, +1), 'z': (5, +1)})
138 mmt2 = (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
139 {'x': (1, +1), 'y': (4, +1), 'z': (6, +1)})
140 mmt3 = (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
141 {'x': (2, +1), 'z': (7, +1)})
142 mmt4 = (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
143 {'x': (8, +1), 'z': (9, +1)})
144 mmt0 = (MomentTensor(m=symmat6(1, 1, 1, 0, 0, 0)),
145 {'x': (0, +1), 'z': (1, +1)})
147 component_scheme = self.store.config.component_scheme
149 if component_scheme == 'elastic8':
150 gfmapping = [mmt1, mmt2, mmt3]
152 elif component_scheme == 'elastic10':
153 gfmapping = [mmt1, mmt2, mmt3, mmt4]
155 elif component_scheme == 'elastic2':
156 gfmapping = [mmt0]
158 elif component_scheme == 'elastic5':
159 gfmapping = [
160 ((1., 1., 0.), {'x': (1, +1), 'z': (4, +1), 'y': (2, +1)}),
161 ((0., 0., 1.), {'x': (0, +1), 'z': (3, +1)})]
162 elif component_scheme == 'elastic18':
163 gfmapping = []
164 for im in range(6):
165 m6 = [0.0] * 6
166 m6[im] = 1.0
168 gfmapping.append(
169 (MomentTensor(m=symmat6(*m6)),
170 {'x': (im, +1),
171 'y': (6+im, +1),
172 'z': (12+im, +1)}))
174 else:
175 raise gf.UnavailableScheme(
176 'fomosto backend "ahfullgreen" cannot handle component scheme '
177 '"%s"' % component_scheme)
179 for source_mech, gfmap in gfmapping:
181 if component_scheme != 'elastic18':
182 norths = xs
183 easts = num.zeros_like(xs)
184 else:
185 receiver = self.store.config.receiver
186 data = []
187 for x in xs:
188 source = clone(self.store.config.source_origin)
189 source.north_shift = x
190 source.east_shift = sy
191 source.depth = sz
192 data.append(receiver.offset_to(source))
194 norths, easts = -num.array(data).T
196 rawtraces = make_traces(
197 self.store.config.earthmodel_1d.require_homogeneous(),
198 source_mech, 1.0/self.store.config.sample_rate,
199 norths, easts, sz, rz, self.ahfullgreen_config.stf)
201 interrupted = []
203 def signal_handler(signum, frame):
204 interrupted.append(True)
206 original = signal.signal(signal.SIGINT, signal_handler)
207 self.store.lock()
208 duplicate_inserts = 0
209 try:
210 for itr, tr in enumerate(rawtraces):
211 if tr.channel not in gfmap:
212 logger.debug('%s not in gfmap' % tr.channel)
213 continue
215 x = xs[tr.meta['isource']]
216 if x > firstx + (nx-1)*dx:
217 logger.error('x out of range')
218 continue
220 ig, factor = gfmap[tr.channel]
222 args = {
223 'A': (sz, x, ig),
224 'B': (rz, sz, x, ig),
225 'C': (sz, sy, x, ig)}[store_type]
227 tr = tr.snap()
229 gf_tr = gf.store.GFTrace.from_trace(tr)
230 gf_tr.data *= factor
232 try:
233 self.store.put(args, gf_tr)
234 except gf.store.DuplicateInsert:
235 duplicate_inserts += 1
237 finally:
238 if duplicate_inserts:
239 logger.warning('%i insertions skipped (duplicates)' %
240 duplicate_inserts)
242 self.store.unlock()
243 signal.signal(signal.SIGINT, original)
245 if interrupted:
246 raise KeyboardInterrupt()
248 logger.info('Done with block %i / %i' %
249 (index+1, self.nblocks))
252def init(store_dir, variant, config_params=None):
253 assert variant is None
255 modelling_code_id = 'ahfullgreen'
257 store_id = os.path.basename(os.path.realpath(store_dir))
259 ahfullgreen = AhfullgreenConfig()
261 d = dict(
262 id=store_id,
263 ncomponents=10,
264 sample_rate=20.,
265 receiver_depth=0*km,
266 source_depth_min=1*km,
267 source_depth_max=10*km,
268 source_depth_delta=1*km,
269 distance_min=1*km,
270 distance_max=20*km,
271 distance_delta=1*km,
272 earthmodel_1d=example_model(),
273 modelling_code_id=modelling_code_id,
274 tabulated_phases=[
275 gf.meta.TPDef(
276 id='anyP',
277 definition='P,p,\\P,\\p'),
278 gf.meta.TPDef(
279 id='anyS',
280 definition='S,s,\\S,\\s')])
282 if config_params is not None:
283 d.update(config_params)
285 config = gf.meta.ConfigTypeA(**d)
286 config.validate()
287 return gf.store.Store.create_editables(
288 store_dir, config=config, extra={'ahfullgreen': ahfullgreen})
291def build(store_dir, force=False, nworkers=None, continue_=False, step=None,
292 iblock=None):
294 return AhfullGFBuilder.build(
295 store_dir, force=force, nworkers=nworkers, continue_=continue_,
296 step=step, iblock=iblock)