1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5from __future__ import absolute_import, division 

6 

7import numpy as num 

8import logging 

9import os 

10import math 

11import signal 

12 

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 

18 

19km = 1000. 

20 

21guts_prefix = 'pf' 

22 

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

24 

25# how to call the programs 

26program_bins = { 

27 'ahfullgreen': 'ahfullgreen', 

28} 

29 

30components = 'r t z'.split() 

31 

32 

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 

40 

41 

42class AhfullgreenError(gf.store.StoreError): 

43 pass 

44 

45 

46class AhfullgreenConfig(Object): 

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

48 

49 

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

51 source_depth, receiver_depth, stf): 

52 

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 

59 

60 npad = 120 

61 

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) 

66 

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

70 

71 outx = num.zeros(ns) 

72 outy = num.zeros(ns) 

73 outz = num.zeros(ns) 

74 

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

76 

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) 

82 

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

90 

91 traces.append(tr) 

92 

93 return traces 

94 

95 

96class AhfullGFBuilder(gf.builder.Builder): 

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

98 

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

100 

101 if block_size is None: 

102 block_size = (1, 1, 2000) 

103 

104 if len(self.store.config.ns) == 2: 

105 block_size = block_size[1:] 

106 

107 gf.builder.Builder.__init__( 

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

109 

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

111 

112 def cleanup(self): 

113 self.store.close() 

114 

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) 

119 

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) 

124 

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

126 (index+1, self.nblocks)) 

127 

128 dx = self.gf_config.distance_delta 

129 

130 distances = num.linspace(firstx, firstx + (nx-1)*dx, nx).tolist() 

131 

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

142 

143 component_scheme = self.store.config.component_scheme 

144 

145 if component_scheme == 'elastic8': 

146 gfmapping = [mmt1, mmt2, mmt3] 

147 

148 elif component_scheme == 'elastic10': 

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

150 

151 elif component_scheme == 'elastic2': 

152 gfmapping = [mmt0] 

153 

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

158 

159 else: 

160 raise gf.UnavailableScheme( 

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

162 '"%s"' % component_scheme) 

163 

164 for source_mech, gfmap in gfmapping: 

165 

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) 

171 

172 interrupted = [] 

173 

174 def signal_handler(signum, frame): 

175 interrupted.append(True) 

176 

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 

185 

186 x = tr.meta['distance'] 

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

188 logger.error("x out of range") 

189 continue 

190 

191 ig, factor = gfmap[tr.channel] 

192 

193 if len(self.store.config.ns) == 2: 

194 args = (sz, x, ig) 

195 else: 

196 args = (rz, sz, x, ig) 

197 

198 tr = tr.snap() 

199 

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

201 gf_tr.data *= factor 

202 

203 try: 

204 self.store.put(args, gf_tr) 

205 except gf.store.DuplicateInsert: 

206 duplicate_inserts += 1 

207 

208 finally: 

209 if duplicate_inserts: 

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

211 duplicate_inserts) 

212 

213 self.store.unlock() 

214 signal.signal(signal.SIGINT, original) 

215 

216 if interrupted: 

217 raise KeyboardInterrupt() 

218 

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

220 (index+1, self.nblocks)) 

221 

222 

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

224 assert variant is None 

225 

226 modelling_code_id = 'ahfullgreen' 

227 

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

229 

230 ahfullgreen = AhfullgreenConfig() 

231 

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

252 

253 if config_params is not None: 

254 d.update(config_params) 

255 

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

257 config.validate() 

258 return gf.store.Store.create_editables( 

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

260 

261 

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

263 iblock=None): 

264 

265 return AhfullGFBuilder.build( 

266 store_dir, force=force, nworkers=nworkers, continue_=continue_, 

267 step=step, iblock=iblock)