Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/fomosto/ahfullgreen.py: 97%

1# http://pyrocko.org - GPLv3

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]

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)