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

1# http://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

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

6import logging

7import os

8import shutil

9from tempfile import mkdtemp

10from subprocess import Popen, PIPE

11from os.path import join as pjoin

13import numpy as num

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

16from pyrocko.guts_array import Array

17from pyrocko import trace, util, gf

19guts_prefix = 'pf'

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

23# how to call the programs

24program_bins = {

25 'poel': 'poel',

26}

29def have_backend():

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

31 try:

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

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

35 except OSError:

36 return False

38 return True

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

44def str_float_vals(vals):

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

48def str_int_vals(vals):

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

52def str_str_vals(vals):

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

56def str_complex_vals(vals):

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

60class PoelSourceFunction(Object):

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

63 def string_for_config(self):

64 return '\n'.join([

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

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

69class PoelModel(Object):

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

72 def string_for_config(self):

73 srows = []

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

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

77 return '\n'.join(srows)

79 def get_nlines(self):

80 return self.data.shape[0]

83class PoelConfig(Object):

84 s_radius = Float.T(default=0.0)

85 s_type = Int.T(default=0)

86 source_function_p = Float.T(default=1.0)

87 source_function_i = PoelSourceFunction.T(

88 default=PoelSourceFunction.D(

89 data=num.array(

90 [[0., 0.], [10., 1.]], dtype=float)))

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

93 accuracy = Float.T(default=0.025)

94 isurfcon = Int.T(default=1)

95 model = PoelModel.T(

96 default=PoelModel.D(

97 data=num.array(

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

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)