1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5# -*- coding: utf-8 -*- 

6from __future__ import absolute_import, division 

7 

8import logging 

9import os 

10import shutil 

11from tempfile import mkdtemp 

12from subprocess import Popen, PIPE 

13from os.path import join as pjoin 

14 

15import numpy as num 

16 

17from pyrocko.guts import Object, Float, Int, List, String 

18from pyrocko.guts_array import Array 

19from pyrocko import trace, util, gf 

20 

21guts_prefix = 'pf' 

22 

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

24 

25# how to call the programs 

26program_bins = { 

27 'poel': 'poel', 

28} 

29 

30 

31def have_backend(): 

32 for cmd in [[exe] for exe in program_bins.values()]: 

33 try: 

34 p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE) 

35 (stdout, stderr) = p.communicate() 

36 

37 except OSError: 

38 return False 

39 

40 return True 

41 

42 

43poel_components = 'uz ur ezz err ett ezr tlt pp dvz dvr'.split() 

44 

45 

46def str_float_vals(vals): 

47 return ' '.join(['%e' % val for val in vals]) 

48 

49 

50def str_int_vals(vals): 

51 return ' '.join(['%i' % val for val in vals]) 

52 

53 

54def str_str_vals(vals): 

55 return ' '.join(["'%s'" % val for val in vals]) 

56 

57 

58def str_complex_vals(vals): 

59 return ', '.join(['(%e, %e)' % (val.real, val.imag) for val in vals]) 

60 

61 

62class PoelSourceFunction(Object): 

63 data = Array.T(shape=(None, 2), dtype=float) 

64 

65 def string_for_config(self): 

66 return '\n'.join([ 

67 '%i %s' % (i+1, str_float_vals(row)) 

68 for (i, row) in enumerate(self.data)]) 

69 

70 

71class PoelModel(Object): 

72 data = Array.T(shape=(None, 6), dtype=float) 

73 

74 def string_for_config(self): 

75 srows = [] 

76 for i, row in enumerate(self.data): 

77 srows.append('%i %s' % (i+1, str_float_vals(row))) 

78 

79 return '\n'.join(srows) 

80 

81 def get_nlines(self): 

82 return self.data.shape[0] 

83 

84 

85class PoelConfig(Object): 

86 s_radius = Float.T(default=0.0) 

87 s_type = Int.T(default=0) 

88 source_function_p = Float.T(default=1.0) 

89 source_function_i = PoelSourceFunction.T( 

90 default=PoelSourceFunction.D( 

91 data=num.array([[0., 0.], [10., 1.]], dtype=float))) 

92 

93 t_window = Float.T(default=500.) 

94 accuracy = Float.T(default=0.025) 

95 isurfcon = Int.T(default=1) 

96 model = PoelModel.T( 

97 default=PoelModel(data=num.array([[ 

98 0.00, 0.4E+09, 0.2, 0.4, 0.75, 5.00]], dtype=float))) 

99 

100 def items(self): 

101 return dict(self.T.inamevals(self)) 

102 

103 

104class PoelConfigFull(PoelConfig): 

105 s_start_depth = Float.T(default=25.0) 

106 s_end_depth = Float.T(default=25.0) 

107 

108 sw_equidistant_z = Int.T(default=1) 

109 no_depths = Int.T(default=10) 

110 depths = Array.T( 

111 shape=(None,), 

112 dtype=float, 

113 default=num.array([10.0, 100.0], dtype=float)) 

114 sw_equidistant_x = Int.T(default=1) 

115 no_distances = Int.T(default=10) 

116 distances = Array.T( 

117 shape=(None,), 

118 dtype=float, 

119 default=num.array([10., 100.])) 

120 

121 no_t_samples = Int.T(default=51) 

122 t_files = List.T(String.T(), default=[x+'.t' for x in poel_components]) 

123 sw_t_files = List.T(Int.T(), default=[1 for x in poel_components]) 

124 

125 def get_output_filenames(self, rundir): 

126 return [pjoin(rundir, fn) for fn in self.t_files] 

127 

128 def string_for_config(self): 

129 

130 d = self.__dict__.copy() 

131 

132 if not self.sw_equidistant_x: 

133 d['no_distances'] = len(self.distances) 

134 d['str_distances'] = str_float_vals(self.distances) 

135 

136 if not self.sw_equidistant_z: 

137 d['no_depths'] = len(self.depths) 

138 d['str_depths'] = str_float_vals(self.depths) 

139 

140 d['sw_t_files_1_2'] = ' '.join( 

141 ['%i' % i for i in self.sw_t_files[0:2]]) 

142 d['t_files_1_2'] = ' '.join( 

143 ["'%s'" % s for s in self.t_files[0:2]]) 

144 d['sw_t_files_3_7'] = ' '.join( 

145 ['%i' % i for i in self.sw_t_files[2:7]]) 

146 d['t_files_3_7'] = ' '.join( 

147 ["'%s'" % s for s in self.t_files[2:7]]) 

148 d['sw_t_files_8_10'] = ' '.join( 

149 ['%i' % i for i in self.sw_t_files[7:10]]) 

150 d['t_files_8_10'] = ' '.join( 

151 ["'%s'" % s for s in self.t_files[7:10]]) 

152 

153 d['no_model_lines'] = self.model.get_nlines() 

154 

155 if self.s_type == 0: 

156 d['source_function'] = str(self.source_function_p) 

157 elif self.s_type == 1: 

158 d['source_function'] = self.source_function_i.string_for_config() 

159 

160 d['model'] = self.model.string_for_config() 

161 

162 template = ''' 

163# This is the input file of FORTRAN77 program "poel06" for modeling 

164# coupled deformation-diffusion processes based on a multi-layered (half- 

165# or full-space) poroelastic media induced by an injection (pump) of 

166# from a borehole or by a (point) reservoir loading. 

167# 

168# by R. Wang, 

169# GeoForschungsZentrum Potsdam 

170# e-mail: wang@gfz-potsdam.de 

171# phone 0049 331 2881209 

172# fax 0049 331 2881204 

173# 

174# Last modified: Potsdam, July, 2012 

175# 

176############################################################## 

177## ## 

178## Cylindrical coordinates (Z positive downwards!) are used ## 

179## If not others specified, SI Unit System is used overall! ## 

180## ## 

181## Tilt is positive when the upper end of a borehole tilt- ## 

182## meter body moves away from the pumping well. ## 

183## ## 

184############################################################## 

185# 

186############################################################################### 

187# 

188# SOURCE PARAMETERS A: SOURCE GEOMETRY 

189# ==================================== 

190# 1. source top and bottom depth [m] 

191# Note: top depth < bottom depth for a vertical line source 

192# top depth = bottom depth for a vertical point source 

193# 

194# ! whole source screen should be within a homogeneous layer, and ! 

195# ! both top and bottom should not coincide with any interface of ! 

196# ! the model used (see below) ! 

197# 

198# 2. source radius (> 0) [m] 

199# Note: source radius > 0 for a horizontal disk source 

200# source radius = 0 for a horizontal point source 

201#------------------------------------------------------------------------------ 

202 %(s_start_depth)g %(s_end_depth)g |dble: s_top_depth, s_bottom_de 

203 %(s_radius)g |dble: s_radius; 

204#------------------------------------------------------------------------------ 

205# 

206# SOURCE PARAMETERS B: SOURCE TYPE 

207# ================================ 

208# 1. selection of source type: 

209# 0 = initial excess pore pressure within the source volume 

210# (initial value problem) 

211# 1 = injection within the source volume 

212# (boundary value problem) 

213#------------------------------------------------------------------------------ 

214 %(s_type)i |int: sw_source_type; 

215#------------------------------------------------------------------------------ 

216 %(source_function)s 

217############################################################################### 

218# 

219# RECEIVER PARAMETERS A: RECEIVER DEPTH SAMPLING 

220# ============================================== 

221# 1. switch for equidistant steping (1/0 = yes/no) 

222# 2. number of receiver depth samples (<= nzrmax defined in "peglobal.h") 

223# 3. if equidistant, start depth [m], end depth [m]; else list of depths 

224# (all >= 0 and ordered from small to large!) 

225#------------------------------------------------------------------------------ 

226 %(sw_equidistant_z)i |int: sw_receiver_depth_samplin 

227 %(no_depths)i |int: no_depths; 

228 %(str_depths)s |dble: zr_1,zr_n; or zr_1,zr_2, 

229#------------------------------------------------------------------------------ 

230# 

231# RECEIVER PARAMETERS B: RECEIVER DISTANCE SAMPLING 

232# ================================================= 

233# 1. switch for equidistant steping (1/0 = yes/no) 

234# 2. number of receiver distance samples (<= nrmax defined in "peglobal.h") 

235# 3. if equidistant, start distance [m], end distance [m]; else list of 

236# distances (all >= 0 and ordered from small to large!) 

237#------------------------------------------------------------------------------ 

238 %(sw_equidistant_x)i |int: sw_equidistant; 

239 %(no_distances)i |int: no_distances; 

240 %(str_distances)s |dble: d_1,d_n; or d_1,d_2, ... 

241#------------------------------------------------------------------------------ 

242# 

243# RECEIVER PARAMETERS C: Time SAMPLING 

244# ==================================== 

245# 1. time window [s] 

246# 2. number of time samples 

247# Note: the caracteristic diffusion time = 

248# max_receiver_distance^2 / diffusivity_of_source_layer 

249#------------------------------------------------------------------------------ 

250 %(t_window)s |dble: time_window; 

251 %(no_t_samples)i |int: no_time_samples; 

252#------------------------------------------------------------------------------ 

253# 

254# WAVENUMBER INTEGRATION PARAMETERS 

255# ================================= 

256# 1. relative accuracy (0.01 for 1%% error) for numerical wavenumber integratio 

257#------------------------------------------------------------------------------ 

258 %(accuracy)s |dble: accuracy; 

259#------------------------------------------------------------------------------ 

260############################################################################### 

261# 

262# OUTPUTS A: DISPLACEMENT TIME SERIES 

263# =================================== 

264# 1. select the 2 displacement time series (1/0 = yes/no) 

265# Note Ut = 0 

266# 2. file names of these 2 time series 

267#------------------------------------------------------------------------------ 

268 %(sw_t_files_1_2)s |int: sw_t_files(1-2); 

269 %(t_files_1_2)s |char: t_files(1-2); 

270#------------------------------------------------------------------------------ 

271# 

272# OUTPUTS B: STRAIN TENSOR & TILT TIME SERIES 

273# =========================================== 

274# 1. select strain time series (1/0 = yes/no): Ezz, Err, Ett, Ezr (4 tensor 

275# components) and Tlt (= -dur/dz, the radial component of the vertical tilt) 

276# Note Ezt, Ert and Tlt (tangential tilt) = 0 

277# 2. file names of these 5 time series 

278#------------------------------------------------------------------------------ 

279 %(sw_t_files_3_7)s |int: sw_t_files(3-7); 

280 %(t_files_3_7)s |char: t_files(3-7); 

281#------------------------------------------------------------------------------ 

282# 

283# OUTPUTS C: PORE PRESSURE & DARCY VELOCITY TIME SERIES 

284# ===================================================== 

285# 1. select pore pressure and Darcy velocity time series (1/0 = yes/no): 

286# Pp (excess pore pressure), Dvz, Dvr (2 Darcy velocity components) 

287# Note Dvt = 0 

288# 2. file names of these 3 time series 

289#------------------------------------------------------------------------------ 

290 %(sw_t_files_8_10)s |int: sw_t_files(8-10); 

291 %(t_files_8_10)s |char: t_files(8-10); 

292#------------------------------------------------------------------------------ 

293# 

294# OUTPUTS D: SNAPSHOTS OF ALL OBSERVABLES 

295# ======================================= 

296# 1. number of snapshots 

297# 2. time[s] (within the time window, see above) and output filename of 

298# the 1. snapshot 

299# 3. ... 

300#------------------------------------------------------------------------------ 

301 1 |int: no_sn; 

302 %(t_window)s 'snapshot.dat' |dable: sn_time(i),sn_file(i) 

303############################################################################### 

304# 

305# GLOBAL MODEL PARAMETERS 

306# ======================= 

307# 1. switch for surface conditions: 

308# 0 = without free surface (whole space), 

309# 1 = unconfined free surface (p = 0), 

310# 2 = confined free surface (dp/dz = 0). 

311# 2. number of data lines of the layered model (<= lmax as defined in 

312# "peglobal.h") (see Note below) 

313#------------------------------------------------------------------------------ 

314 %(isurfcon)i |int: isurfcon 

315 %(no_model_lines)i |int: no_model_lines; 

316#------------------------------------------------------------------------------ 

317# 

318# MULTILAYERED MODEL PARAMETERS 

319# ============================= 

320# 

321# Note: mu = shear modulus 

322# nu = Poisson ratio under drained condition 

323# nu_u = Poisson ratio under undrained condition (nu_u > nu) 

324# B = Skempton parameter (the change in pore pressure per unit change 

325# in confining pressure under undrained condition) 

326# D = hydraulic diffusivity 

327# 

328# no depth[m] mu[Pa] nu nu_u B D[m^2/s] Explanations 

329#------------------------------------------------------------------------------ 

330 %(model)s 

331###########################end of all inputs################################### 

332 

333Note for the model input format and the step-function approximation for model 

334parameters varying linearly with depth: 

335 

336The surface and the upper boundary of the lowest half-space as well as the 

337interfaces at which the poroelastic parameters are continuous, are all defined 

338by a single data line; All other interfaces, at which the poroelastic 

339parameters are discontinuous, are all defined by two data lines (upper-side and 

340lower-side values). This input format would also be needed for a graphic plot 

341of the layered model. Layers which have different parameter values at top and 

342bottom, will be treated as layers with a constant gradient, and will be 

343discretised to a number of homogeneous sublayers. Errors due to the 

344discretisation are limited within about 5%% (changeable, see peglobal.h). 

345'''.lstrip() 

346 

347 return template % d 

348 

349 

350class PoelError(Exception): 

351 pass 

352 

353 

354class PoelRunner(object): 

355 

356 def __init__(self, tmp=None): 

357 

358 self.tempdir = mkdtemp(prefix='poelrun', dir=tmp) 

359 self.program = program_bins['poel'] 

360 self.config = None 

361 

362 def run(self, config): 

363 self.config = config 

364 

365 input_fn = pjoin(self.tempdir, 'input') 

366 

367 f = open(input_fn, 'w') 

368 poel_input = config.string_for_config() 

369 

370 logger.debug( 

371 '===== begin poel input =====\n%s===== end poel input =====' 

372 % poel_input) 

373 

374 f.write(poel_input) 

375 f.close() 

376 program = self.program 

377 

378 old_wd = os.getcwd() 

379 

380 os.chdir(self.tempdir) 

381 

382 try: 

383 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE) 

384 except OSError: 

385 os.chdir(old_wd) 

386 raise PoelError( 

387 '''could not start poel executable: "%s" 

388Available fomosto backends and download links to the modelling codes are listed 

389on 

390 

391 https://pyrocko.org/docs/current/apps/fomosto/backends.html 

392 

393''' % program) 

394 

395 (poel_output, poel_error) = proc.communicate('input\n') 

396 

397 logger.debug( 

398 '===== begin poel output =====\n%s===== end poel output =====' 

399 % poel_output) 

400 if poel_error: 

401 logger.error( 

402 '===== begin poel error =====\n%s===== end poel error =====' 

403 % poel_error) 

404 

405 errmess = [] 

406 if proc.returncode != 0: 

407 errmess.append( 

408 'poel had a non-zero exit state: %i' % proc.returncode) 

409 if poel_error: 

410 errmess.append('poel emitted something via stderr') 

411 if poel_output.lower().find('error') != -1: 

412 errmess.append("the string 'error' appeared in poel output") 

413 

414 if errmess: 

415 os.chdir(old_wd) 

416 raise PoelError( 

417 '''===== begin poel input =====\n%s===== end poel input ===== 

418===== begin poel output =====\n%s===== end poel output ===== 

419===== begin poel error =====\n%s===== end poel error ===== 

420%s 

421poel has been invoked as "%s"''' % ( 

422 poel_input, 

423 poel_output, 

424 poel_error, 

425 '\n'.join(errmess), 

426 program)) 

427 

428 self.poel_output = poel_output 

429 self.poel_error = poel_error 

430 

431 os.chdir(old_wd) 

432 

433 def get_traces(self): 

434 

435 if self.config.sw_equidistant_x == 1: 

436 nx = self.config.no_distances 

437 xmin, xmax = self.config.distances 

438 if nx > 1: 

439 dx = (xmax-xmin)/(nx-1) 

440 else: 

441 dx = 1.0 

442 distances = [xmin + ix*dx for ix in range(nx)] 

443 else: 

444 distances = self.config.distances 

445 

446 if self.config.sw_equidistant_z == 1: 

447 nrz = self.config.no_depths 

448 rzmin, rzmax = self.config.depths 

449 if nrz > 1: 

450 drz = (rzmax-rzmin)/(nrz-1) 

451 else: 

452 drz = 1.0 

453 rdepths = [rzmin + irz*drz for irz in range(nrz)] 

454 else: 

455 rdepths = self.config.depths 

456 

457 sz = self.config.s_start_depth 

458 

459 fns = self.config.get_output_filenames(self.tempdir) 

460 traces = [] 

461 for comp, fn in zip(poel_components, fns): 

462 data = num.loadtxt(fn, skiprows=1, dtype=float) 

463 nsamples, ntraces = data.shape 

464 ntraces -= 1 

465 tmin = data[0, 0] 

466 deltat = (data[-1, 0] - data[0, 0])/(nsamples-1) 

467 for itrace in range(ntraces): 

468 x = distances[itrace % len(distances)] 

469 rz = rdepths[itrace // len(distances)] 

470 tr = trace.Trace( 

471 '', '%i' % itrace, 'c', comp, 

472 tmin=tmin, deltat=deltat, ydata=data[:, itrace+1], 

473 meta={'itrace': itrace, 'x': x, 'rz': rz, 'sz': sz}) 

474 

475 traces.append(tr) 

476 

477 return traces 

478 

479 def __del__(self): 

480 shutil.rmtree(self.tempdir) 

481 

482 

483class PoelGFBuilder(gf.builder.Builder): 

484 def __init__(self, store_dir, step, shared, block_size=None, tmp=None, 

485 force=False): 

486 

487 if block_size is None: 

488 block_size = (51, 1, 51) 

489 

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

491 

492 gf.builder.Builder.__init__( 

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

494 

495 self.poel_config = self.store.get_extra('poel') 

496 

497 self.tmp = tmp 

498 if self.tmp is not None: 

499 util.ensuredir(self.tmp) 

500 

501 def cleanup(self): 

502 self.store.close() 

503 

504 def work_block(self, index): 

505 logger.info('Starting block %i / %i' % (index+1, self.nblocks)) 

506 

507 runner = PoelRunner(tmp=self.tmp) 

508 

509 conf = PoelConfigFull(**self.poel_config.items()) 

510 

511 (firstrz, sz, firstx), (lastrz, sz, lastx), (nrz, _, nx) = \ 

512 self.get_block_extents(index) 

513 

514 conf.s_start_depth = sz 

515 conf.s_end_depth = sz 

516 conf.sw_equidistant_x = 1 

517 conf.distances = [firstx, lastx] 

518 conf.sw_equidistant_z = 1 

519 conf.no_distances = nx 

520 conf.depths = [firstrz, lastrz] 

521 conf.no_depths = nrz 

522 conf.no_t_samples = int( 

523 round(conf.t_window * self.gf_config.sample_rate)) + 1 

524 

525 runner.run(conf) 

526 

527 comp2ig = dict([(c, ig) for (ig, c) in enumerate(poel_components)]) 

528 

529 rawtraces = runner.get_traces() 

530 

531 self.store.lock() 

532 

533 for tr in rawtraces: 

534 

535 x = tr.meta['x'] 

536 rz = tr.meta['rz'] 

537 sz = tr.meta['sz'] 

538 

539 ig = comp2ig[tr.channel] 

540 

541 gf_tr = gf.store.GFTrace( 

542 tr.get_ydata(), 

543 int(round(tr.tmin / tr.deltat)), 

544 tr.deltat) 

545 

546 self.store.put((rz, sz, x, ig), gf_tr) 

547 

548 self.store.unlock() 

549 

550 logger.info('Done with block %i / %i' % (index+1, self.nblocks)) 

551 

552 

553def init(store_dir, variant): 

554 if variant is not None: 

555 raise gf.store.StoreError('unsupported variant: %s' % variant) 

556 

557 poel = PoelConfig() 

558 

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

560 

561 config = gf.meta.ConfigTypeB( 

562 modelling_code_id='poel', 

563 id=store_id, 

564 ncomponents=10, 

565 component_scheme='poroelastic10', 

566 sample_rate=0.1, 

567 distance_min=10.0, 

568 distance_max=20.0, 

569 distance_delta=5.0, 

570 source_depth_min=10.0, 

571 source_depth_max=20.0, 

572 source_depth_delta=5.0, 

573 receiver_depth_min=10.0, 

574 receiver_depth_max=10.0, 

575 receiver_depth_delta=5.0) 

576 

577 return gf.store.Store.create_editables( 

578 store_dir, 

579 config=config, 

580 extra={'poel': poel}) 

581 

582 

583def build( 

584 store_dir, 

585 force=False, 

586 nworkers=None, 

587 continue_=False, 

588 step=None, 

589 iblock=None): 

590 

591 return PoelGFBuilder.build( 

592 store_dir, 

593 force=force, 

594 nworkers=nworkers, 

595 continue_=continue_, 

596 step=step, 

597 iblock=iblock)