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
8import logging
9import os
10import shutil
11from tempfile import mkdtemp
12from subprocess import Popen, PIPE
13from os.path import join as pjoin
15import numpy as num
17from pyrocko.guts import Object, Float, Int, List, String
18from pyrocko.guts_array import Array
19from pyrocko import trace, util, gf
21guts_prefix = 'pf'
23logger = logging.getLogger('pyrocko.fomosto.poel')
25# how to call the programs
26program_bins = {
27 'poel': 'poel',
28}
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()
37 except OSError:
38 return False
40 return True
43poel_components = 'uz ur ezz err ett ezr tlt pp dvz dvr'.split()
46def str_float_vals(vals):
47 return ' '.join(['%e' % val for val in vals])
50def str_int_vals(vals):
51 return ' '.join(['%i' % val for val in vals])
54def str_str_vals(vals):
55 return ' '.join(["'%s'" % val for val in vals])
58def str_complex_vals(vals):
59 return ', '.join(['(%e, %e)' % (val.real, val.imag) for val in vals])
62class PoelSourceFunction(Object):
63 data = Array.T(shape=(None, 2), dtype=float)
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)])
71class PoelModel(Object):
72 data = Array.T(shape=(None, 6), dtype=float)
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)))
79 return '\n'.join(srows)
81 def get_nlines(self):
82 return self.data.shape[0]
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)))
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)))
100 def items(self):
101 return dict(self.T.inamevals(self))
104class PoelConfigFull(PoelConfig):
105 s_start_depth = Float.T(default=25.0)
106 s_end_depth = Float.T(default=25.0)
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.]))
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])
125 def get_output_filenames(self, rundir):
126 return [pjoin(rundir, fn) for fn in self.t_files]
128 def string_for_config(self):
130 d = self.__dict__.copy()
132 if not self.sw_equidistant_x:
133 d['no_distances'] = len(self.distances)
134 d['str_distances'] = str_float_vals(self.distances)
136 if not self.sw_equidistant_z:
137 d['no_depths'] = len(self.depths)
138 d['str_depths'] = str_float_vals(self.depths)
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]])
153 d['no_model_lines'] = self.model.get_nlines()
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()
160 d['model'] = self.model.string_for_config()
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###################################
333Note for the model input format and the step-function approximation for model
334parameters varying linearly with depth:
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()
347 return template % d
350class PoelError(Exception):
351 pass
354class PoelRunner(object):
356 def __init__(self, tmp=None):
358 self.tempdir = mkdtemp(prefix='poelrun', dir=tmp)
359 self.program = program_bins['poel']
360 self.config = None
362 def run(self, config):
363 self.config = config
365 input_fn = pjoin(self.tempdir, 'input')
367 f = open(input_fn, 'w')
368 poel_input = config.string_for_config()
370 logger.debug(
371 '===== begin poel input =====\n%s===== end poel input ====='
372 % poel_input)
374 f.write(poel_input)
375 f.close()
376 program = self.program
378 old_wd = os.getcwd()
380 os.chdir(self.tempdir)
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
391 https://pyrocko.org/docs/current/apps/fomosto/backends.html
393''' % program)
395 (poel_output, poel_error) = proc.communicate('input\n')
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)
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")
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))
428 self.poel_output = poel_output
429 self.poel_error = poel_error
431 os.chdir(old_wd)
433 def get_traces(self):
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
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
457 sz = self.config.s_start_depth
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})
475 traces.append(tr)
477 return traces
479 def __del__(self):
480 shutil.rmtree(self.tempdir)
483class PoelGFBuilder(gf.builder.Builder):
484 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
485 force=False):
487 if block_size is None:
488 block_size = (51, 1, 51)
490 self.store = gf.store.Store(store_dir, 'w')
492 gf.builder.Builder.__init__(
493 self, self.store.config, step, block_size=block_size, force=force)
495 self.poel_config = self.store.get_extra('poel')
497 self.tmp = tmp
498 if self.tmp is not None:
499 util.ensuredir(self.tmp)
501 def cleanup(self):
502 self.store.close()
504 def work_block(self, index):
505 logger.info('Starting block %i / %i' % (index+1, self.nblocks))
507 runner = PoelRunner(tmp=self.tmp)
509 conf = PoelConfigFull(**self.poel_config.items())
511 (firstrz, sz, firstx), (lastrz, sz, lastx), (nrz, _, nx) = \
512 self.get_block_extents(index)
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
525 runner.run(conf)
527 comp2ig = dict([(c, ig) for (ig, c) in enumerate(poel_components)])
529 rawtraces = runner.get_traces()
531 self.store.lock()
533 for tr in rawtraces:
535 x = tr.meta['x']
536 rz = tr.meta['rz']
537 sz = tr.meta['sz']
539 ig = comp2ig[tr.channel]
541 gf_tr = gf.store.GFTrace(
542 tr.get_ydata(),
543 int(round(tr.tmin / tr.deltat)),
544 tr.deltat)
546 self.store.put((rz, sz, x, ig), gf_tr)
548 self.store.unlock()
550 logger.info('Done with block %i / %i' % (index+1, self.nblocks))
553def init(store_dir, variant):
554 if variant is not None:
555 raise gf.store.StoreError('unsupported variant: %s' % variant)
557 poel = PoelConfig()
559 store_id = os.path.basename(os.path.realpath(store_dir))
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)
577 return gf.store.Store.create_editables(
578 store_dir,
579 config=config,
580 extra={'poel': poel})
583def build(
584 store_dir,
585 force=False,
586 nworkers=None,
587 continue_=False,
588 step=None,
589 iblock=None):
591 return PoelGFBuilder.build(
592 store_dir,
593 force=force,
594 nworkers=nworkers,
595 continue_=continue_,
596 step=step,
597 iblock=iblock)