1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6import numpy as num 

7import logging 

8import os 

9import math 

10import signal 

11 

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 

17 

18km = 1000. 

19 

20guts_prefix = 'pf' 

21 

22logger = logging.getLogger('pyrocko.fomosto.ahfullgreen') 

23 

24# how to call the programs 

25program_bins = { 

26 'ahfullgreen': 'ahfullgreen', 

27} 

28 

29components = 'x y z'.split() 

30 

31 

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 

39 

40 

41class AhfullgreenError(gf.store.StoreError): 

42 pass 

43 

44 

45class AhfullgreenConfig(Object): 

46 stf = AhfullgreenSTF.T(default=AhfullgreenSTFImpulse.D()) 

47 

48 

49def make_traces(material, source_mech, deltat, norths, easts, 

50 source_depth, receiver_depth, stf): 

51 

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 

58 

59 npad = 120 

60 

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) 

65 

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)) 

69 

70 outx = num.zeros(ns) 

71 outy = num.zeros(ns) 

72 outz = num.zeros(ns) 

73 

74 x = (north, east, receiver_depth-source_depth) 

75 

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) 

81 

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)) 

87 

88 traces.append(tr) 

89 

90 return traces 

91 

92 

93class AhfullGFBuilder(gf.builder.Builder): 

94 def __init__(self, store_dir, step, shared, block_size=None, force=False): 

95 

96 self.store = gf.store.Store(store_dir, 'w') 

97 

98 block_size = { 

99 'A': (1, 2000), 

100 'B': (1, 1, 2000), 

101 'C': (1, 1, 2)}[self.store.config.short_type] 

102 

103 gf.builder.Builder.__init__( 

104 self, self.store.config, step, block_size=block_size, force=force) 

105 

106 self.ahfullgreen_config = self.store.get_extra('ahfullgreen') 

107 

108 def cleanup(self): 

109 self.store.close() 

110 

111 def work_block(self, index): 

112 store_type = self.store.config.short_type 

113 

114 if store_type == 'A': 

115 (sz, firstx), (sz, lastx), (ns, nx) = \ 

116 self.get_block_extents(index) 

117 

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 

128 

129 logger.info('Starting block %i / %i' % 

130 (index+1, self.nblocks)) 

131 

132 dx = self.gf_config.deltas[-1] 

133 

134 xs = num.linspace(firstx, firstx + (nx-1)*dx, nx).tolist() 

135 

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)}) 

146 

147 component_scheme = self.store.config.component_scheme 

148 

149 if component_scheme == 'elastic8': 

150 gfmapping = [mmt1, mmt2, mmt3] 

151 

152 elif component_scheme == 'elastic10': 

153 gfmapping = [mmt1, mmt2, mmt3, mmt4] 

154 

155 elif component_scheme == 'elastic2': 

156 gfmapping = [mmt0] 

157 

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 

167 

168 gfmapping.append( 

169 (MomentTensor(m=symmat6(*m6)), 

170 {'x': (im, +1), 

171 'y': (6+im, +1), 

172 'z': (12+im, +1)})) 

173 

174 else: 

175 raise gf.UnavailableScheme( 

176 'fomosto backend "ahfullgreen" cannot handle component scheme ' 

177 '"%s"' % component_scheme) 

178 

179 for source_mech, gfmap in gfmapping: 

180 

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)) 

193 

194 norths, easts = -num.array(data).T 

195 

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) 

200 

201 interrupted = [] 

202 

203 def signal_handler(signum, frame): 

204 interrupted.append(True) 

205 

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 

214 

215 x = xs[tr.meta['isource']] 

216 if x > firstx + (nx-1)*dx: 

217 logger.error('x out of range') 

218 continue 

219 

220 ig, factor = gfmap[tr.channel] 

221 

222 args = { 

223 'A': (sz, x, ig), 

224 'B': (rz, sz, x, ig), 

225 'C': (sz, sy, x, ig)}[store_type] 

226 

227 tr = tr.snap() 

228 

229 gf_tr = gf.store.GFTrace.from_trace(tr) 

230 gf_tr.data *= factor 

231 

232 try: 

233 self.store.put(args, gf_tr) 

234 except gf.store.DuplicateInsert: 

235 duplicate_inserts += 1 

236 

237 finally: 

238 if duplicate_inserts: 

239 logger.warning('%i insertions skipped (duplicates)' % 

240 duplicate_inserts) 

241 

242 self.store.unlock() 

243 signal.signal(signal.SIGINT, original) 

244 

245 if interrupted: 

246 raise KeyboardInterrupt() 

247 

248 logger.info('Done with block %i / %i' % 

249 (index+1, self.nblocks)) 

250 

251 

252def init(store_dir, variant, config_params=None): 

253 assert variant is None 

254 

255 modelling_code_id = 'ahfullgreen' 

256 

257 store_id = os.path.basename(os.path.realpath(store_dir)) 

258 

259 ahfullgreen = AhfullgreenConfig() 

260 

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')]) 

281 

282 if config_params is not None: 

283 d.update(config_params) 

284 

285 config = gf.meta.ConfigTypeA(**d) 

286 config.validate() 

287 return gf.store.Store.create_editables( 

288 store_dir, config=config, extra={'ahfullgreen': ahfullgreen}) 

289 

290 

291def build(store_dir, force=False, nworkers=None, continue_=False, step=None, 

292 iblock=None): 

293 

294 return AhfullGFBuilder.build( 

295 store_dir, force=force, nworkers=nworkers, continue_=continue_, 

296 step=step, iblock=iblock)