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
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 = 'r t 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 i_distance, (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' % i_distance, '', comp,
86 tmin=tmin, ydata=o, deltat=deltat,
87 meta=dict(
88 distance=math.sqrt(north**2 + east**2),
89 azimuth=0.))
91 traces.append(tr)
93 return traces
96class AhfullGFBuilder(gf.builder.Builder):
97 def __init__(self, store_dir, step, shared, block_size=None, force=False):
99 self.store = gf.store.Store(store_dir, 'w')
101 if block_size is None:
102 block_size = (1, 1, 2000)
104 if len(self.store.config.ns) == 2:
105 block_size = block_size[1:]
107 gf.builder.Builder.__init__(
108 self, self.store.config, step, block_size=block_size, force=force)
110 self.ahfullgreen_config = self.store.get_extra('ahfullgreen')
112 def cleanup(self):
113 self.store.close()
115 def work_block(self, index):
116 if len(self.store.config.ns) == 2:
117 (sz, firstx), (sz, lastx), (ns, nx) = \
118 self.get_block_extents(index)
120 rz = self.store.config.receiver_depth
121 else:
122 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
123 self.get_block_extents(index)
125 logger.info('Starting block %i / %i' %
126 (index+1, self.nblocks))
128 dx = self.gf_config.distance_delta
130 distances = num.linspace(firstx, firstx + (nx-1)*dx, nx).tolist()
132 mmt1 = (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
133 {'r': (0, +1), 't': (3, +1), 'z': (5, +1)})
134 mmt2 = (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
135 {'r': (1, +1), 't': (4, +1), 'z': (6, +1)})
136 mmt3 = (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
137 {'r': (2, +1), 'z': (7, +1)})
138 mmt4 = (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
139 {'r': (8, +1), 'z': (9, +1)})
140 mmt0 = (MomentTensor(m=symmat6(1, 1, 1, 0, 0, 0)),
141 {'r': (0, +1), 'z': (1, +1)})
143 component_scheme = self.store.config.component_scheme
145 if component_scheme == 'elastic8':
146 gfmapping = [mmt1, mmt2, mmt3]
148 elif component_scheme == 'elastic10':
149 gfmapping = [mmt1, mmt2, mmt3, mmt4]
151 elif component_scheme == 'elastic2':
152 gfmapping = [mmt0]
154 elif component_scheme == 'elastic5':
155 gfmapping = [
156 ((1., 1., 0.), {'r': (1, +1), 'z': (4, +1), 't': (2, +1)}),
157 ((0., 0., 1.), {'r': (0, +1), 'z': (3, +1)})]
159 else:
160 raise gf.UnavailableScheme(
161 'fomosto backend "ahfullgreen" cannot handle component scheme '
162 '"%s"' % component_scheme)
164 for source_mech, gfmap in gfmapping:
166 rawtraces = make_traces(
167 self.store.config.earthmodel_1d.require_homogeneous(),
168 source_mech, 1.0/self.store.config.sample_rate,
169 distances, num.zeros_like(distances), sz, rz,
170 self.ahfullgreen_config.stf)
172 interrupted = []
174 def signal_handler(signum, frame):
175 interrupted.append(True)
177 original = signal.signal(signal.SIGINT, signal_handler)
178 self.store.lock()
179 duplicate_inserts = 0
180 try:
181 for itr, tr in enumerate(rawtraces):
182 if tr.channel not in gfmap:
183 logger.debug('%s not in gfmap' % tr.channel)
184 continue
186 x = tr.meta['distance']
187 if x > firstx + (nx-1)*dx:
188 logger.error("x out of range")
189 continue
191 ig, factor = gfmap[tr.channel]
193 if len(self.store.config.ns) == 2:
194 args = (sz, x, ig)
195 else:
196 args = (rz, sz, x, ig)
198 tr = tr.snap()
200 gf_tr = gf.store.GFTrace.from_trace(tr)
201 gf_tr.data *= factor
203 try:
204 self.store.put(args, gf_tr)
205 except gf.store.DuplicateInsert:
206 duplicate_inserts += 1
208 finally:
209 if duplicate_inserts:
210 logger.warn('%i insertions skipped (duplicates)' %
211 duplicate_inserts)
213 self.store.unlock()
214 signal.signal(signal.SIGINT, original)
216 if interrupted:
217 raise KeyboardInterrupt()
219 logger.info('Done with block %i / %i' %
220 (index+1, self.nblocks))
223def init(store_dir, variant, config_params=None):
224 assert variant is None
226 modelling_code_id = 'ahfullgreen'
228 store_id = os.path.basename(os.path.realpath(store_dir))
230 ahfullgreen = AhfullgreenConfig()
232 d = dict(
233 id=store_id,
234 ncomponents=10,
235 sample_rate=20.,
236 receiver_depth=0*km,
237 source_depth_min=1*km,
238 source_depth_max=10*km,
239 source_depth_delta=1*km,
240 distance_min=1*km,
241 distance_max=20*km,
242 distance_delta=1*km,
243 earthmodel_1d=example_model(),
244 modelling_code_id=modelling_code_id,
245 tabulated_phases=[
246 gf.meta.TPDef(
247 id='anyP',
248 definition='P,p,\\P,\\p'),
249 gf.meta.TPDef(
250 id='anyS',
251 definition='S,s,\\S,\\s')])
253 if config_params is not None:
254 d.update(config_params)
256 config = gf.meta.ConfigTypeA(**d)
257 config.validate()
258 return gf.store.Store.create_editables(
259 store_dir, config=config, extra={'ahfullgreen': ahfullgreen})
262def build(store_dir, force=False, nworkers=None, continue_=False, step=None,
263 iblock=None):
265 return AhfullGFBuilder.build(
266 store_dir, force=force, nworkers=nworkers, continue_=continue_,
267 step=step, iblock=iblock)