1# http://pyrocko.org - GPLv3
2#
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([[0., 0.], [10., 1.]], dtype=float)))
91 t_window = Float.T(default=500.)
92 accuracy = Float.T(default=0.025)
93 isurfcon = Int.T(default=1)
94 model = PoelModel.T(
95 default=PoelModel(data=num.array([[
96 0.00, 0.4E+09, 0.2, 0.4, 0.75, 5.00]], dtype=float)))
98 def items(self):
99 return dict(self.T.inamevals(self))
102class PoelConfigFull(PoelConfig):
103 s_start_depth = Float.T(default=25.0)
104 s_end_depth = Float.T(default=25.0)
106 sw_equidistant_z = Int.T(default=1)
107 no_depths = Int.T(default=10)
108 depths = Array.T(
109 shape=(None,),
110 dtype=float,
111 default=num.array([10.0, 100.0], dtype=float))
112 sw_equidistant_x = Int.T(default=1)
113 no_distances = Int.T(default=10)
114 distances = Array.T(
115 shape=(None,),
116 dtype=float,
117 default=num.array([10., 100.]))
119 no_t_samples = Int.T(default=51)
120 t_files = List.T(String.T(), default=[x+'.t' for x in poel_components])
121 sw_t_files = List.T(Int.T(), default=[1 for x in poel_components])
123 def get_output_filenames(self, rundir):
124 return [pjoin(rundir, fn) for fn in self.t_files]
126 def string_for_config(self):
128 d = self.__dict__.copy()
130 if not self.sw_equidistant_x:
131 d['no_distances'] = len(self.distances)
132 d['str_distances'] = str_float_vals(self.distances)
134 if not self.sw_equidistant_z:
135 d['no_depths'] = len(self.depths)
136 d['str_depths'] = str_float_vals(self.depths)
138 d['sw_t_files_1_2'] = ' '.join(
139 ['%i' % i for i in self.sw_t_files[0:2]])
140 d['t_files_1_2'] = ' '.join(
141 ["'%s'" % s for s in self.t_files[0:2]])
142 d['sw_t_files_3_7'] = ' '.join(
143 ['%i' % i for i in self.sw_t_files[2:7]])
144 d['t_files_3_7'] = ' '.join(
145 ["'%s'" % s for s in self.t_files[2:7]])
146 d['sw_t_files_8_10'] = ' '.join(
147 ['%i' % i for i in self.sw_t_files[7:10]])
148 d['t_files_8_10'] = ' '.join(
149 ["'%s'" % s for s in self.t_files[7:10]])
151 d['no_model_lines'] = self.model.get_nlines()
153 if self.s_type == 0:
154 d['source_function'] = str(self.source_function_p)
155 elif self.s_type == 1:
156 d['source_function'] = self.source_function_i.string_for_config()
158 d['model'] = self.model.string_for_config()
160 template = '''
161# This is the input file of FORTRAN77 program "poel06" for modeling
162# coupled deformation-diffusion processes based on a multi-layered (half-
163# or full-space) poroelastic media induced by an injection (pump) of
164# from a borehole or by a (point) reservoir loading.
165#
166# by R. Wang,
167# GeoForschungsZentrum Potsdam
168# e-mail: wang@gfz-potsdam.de
169# phone 0049 331 2881209
170# fax 0049 331 2881204
171#
172# Last modified: Potsdam, July, 2012
173#
174##############################################################
175## ##
176## Cylindrical coordinates (Z positive downwards!) are used ##
177## If not others specified, SI Unit System is used overall! ##
178## ##
179## Tilt is positive when the upper end of a borehole tilt- ##
180## meter body moves away from the pumping well. ##
181## ##
182##############################################################
183#
184###############################################################################
185#
186# SOURCE PARAMETERS A: SOURCE GEOMETRY
187# ====================================
188# 1. source top and bottom depth [m]
189# Note: top depth < bottom depth for a vertical line source
190# top depth = bottom depth for a vertical point source
191#
192# ! whole source screen should be within a homogeneous layer, and !
193# ! both top and bottom should not coincide with any interface of !
194# ! the model used (see below) !
195#
196# 2. source radius (> 0) [m]
197# Note: source radius > 0 for a horizontal disk source
198# source radius = 0 for a horizontal point source
199#------------------------------------------------------------------------------
200 %(s_start_depth)g %(s_end_depth)g |dble: s_top_depth, s_bottom_de
201 %(s_radius)g |dble: s_radius;
202#------------------------------------------------------------------------------
203#
204# SOURCE PARAMETERS B: SOURCE TYPE
205# ================================
206# 1. selection of source type:
207# 0 = initial excess pore pressure within the source volume
208# (initial value problem)
209# 1 = injection within the source volume
210# (boundary value problem)
211#------------------------------------------------------------------------------
212 %(s_type)i |int: sw_source_type;
213#------------------------------------------------------------------------------
214 %(source_function)s
215###############################################################################
216#
217# RECEIVER PARAMETERS A: RECEIVER DEPTH SAMPLING
218# ==============================================
219# 1. switch for equidistant steping (1/0 = yes/no)
220# 2. number of receiver depth samples (<= nzrmax defined in "peglobal.h")
221# 3. if equidistant, start depth [m], end depth [m]; else list of depths
222# (all >= 0 and ordered from small to large!)
223#------------------------------------------------------------------------------
224 %(sw_equidistant_z)i |int: sw_receiver_depth_samplin
225 %(no_depths)i |int: no_depths;
226 %(str_depths)s |dble: zr_1,zr_n; or zr_1,zr_2,
227#------------------------------------------------------------------------------
228#
229# RECEIVER PARAMETERS B: RECEIVER DISTANCE SAMPLING
230# =================================================
231# 1. switch for equidistant steping (1/0 = yes/no)
232# 2. number of receiver distance samples (<= nrmax defined in "peglobal.h")
233# 3. if equidistant, start distance [m], end distance [m]; else list of
234# distances (all >= 0 and ordered from small to large!)
235#------------------------------------------------------------------------------
236 %(sw_equidistant_x)i |int: sw_equidistant;
237 %(no_distances)i |int: no_distances;
238 %(str_distances)s |dble: d_1,d_n; or d_1,d_2, ...
239#------------------------------------------------------------------------------
240#
241# RECEIVER PARAMETERS C: Time SAMPLING
242# ====================================
243# 1. time window [s]
244# 2. number of time samples
245# Note: the caracteristic diffusion time =
246# max_receiver_distance^2 / diffusivity_of_source_layer
247#------------------------------------------------------------------------------
248 %(t_window)s |dble: time_window;
249 %(no_t_samples)i |int: no_time_samples;
250#------------------------------------------------------------------------------
251#
252# WAVENUMBER INTEGRATION PARAMETERS
253# =================================
254# 1. relative accuracy (0.01 for 1%% error) for numerical wavenumber integratio
255#------------------------------------------------------------------------------
256 %(accuracy)s |dble: accuracy;
257#------------------------------------------------------------------------------
258###############################################################################
259#
260# OUTPUTS A: DISPLACEMENT TIME SERIES
261# ===================================
262# 1. select the 2 displacement time series (1/0 = yes/no)
263# Note Ut = 0
264# 2. file names of these 2 time series
265#------------------------------------------------------------------------------
266 %(sw_t_files_1_2)s |int: sw_t_files(1-2);
267 %(t_files_1_2)s |char: t_files(1-2);
268#------------------------------------------------------------------------------
269#
270# OUTPUTS B: STRAIN TENSOR & TILT TIME SERIES
271# ===========================================
272# 1. select strain time series (1/0 = yes/no): Ezz, Err, Ett, Ezr (4 tensor
273# components) and Tlt (= -dur/dz, the radial component of the vertical tilt)
274# Note Ezt, Ert and Tlt (tangential tilt) = 0
275# 2. file names of these 5 time series
276#------------------------------------------------------------------------------
277 %(sw_t_files_3_7)s |int: sw_t_files(3-7);
278 %(t_files_3_7)s |char: t_files(3-7);
279#------------------------------------------------------------------------------
280#
281# OUTPUTS C: PORE PRESSURE & DARCY VELOCITY TIME SERIES
282# =====================================================
283# 1. select pore pressure and Darcy velocity time series (1/0 = yes/no):
284# Pp (excess pore pressure), Dvz, Dvr (2 Darcy velocity components)
285# Note Dvt = 0
286# 2. file names of these 3 time series
287#------------------------------------------------------------------------------
288 %(sw_t_files_8_10)s |int: sw_t_files(8-10);
289 %(t_files_8_10)s |char: t_files(8-10);
290#------------------------------------------------------------------------------
291#
292# OUTPUTS D: SNAPSHOTS OF ALL OBSERVABLES
293# =======================================
294# 1. number of snapshots
295# 2. time[s] (within the time window, see above) and output filename of
296# the 1. snapshot
297# 3. ...
298#------------------------------------------------------------------------------
299 1 |int: no_sn;
300 %(t_window)s 'snapshot.dat' |dable: sn_time(i),sn_file(i)
301###############################################################################
302#
303# GLOBAL MODEL PARAMETERS
304# =======================
305# 1. switch for surface conditions:
306# 0 = without free surface (whole space),
307# 1 = unconfined free surface (p = 0),
308# 2 = confined free surface (dp/dz = 0).
309# 2. number of data lines of the layered model (<= lmax as defined in
310# "peglobal.h") (see Note below)
311#------------------------------------------------------------------------------
312 %(isurfcon)i |int: isurfcon
313 %(no_model_lines)i |int: no_model_lines;
314#------------------------------------------------------------------------------
315#
316# MULTILAYERED MODEL PARAMETERS
317# =============================
318#
319# Note: mu = shear modulus
320# nu = Poisson ratio under drained condition
321# nu_u = Poisson ratio under undrained condition (nu_u > nu)
322# B = Skempton parameter (the change in pore pressure per unit change
323# in confining pressure under undrained condition)
324# D = hydraulic diffusivity
325#
326# no depth[m] mu[Pa] nu nu_u B D[m^2/s] Explanations
327#------------------------------------------------------------------------------
328 %(model)s
329###########################end of all inputs###################################
331Note for the model input format and the step-function approximation for model
332parameters varying linearly with depth:
334The surface and the upper boundary of the lowest half-space as well as the
335interfaces at which the poroelastic parameters are continuous, are all defined
336by a single data line; All other interfaces, at which the poroelastic
337parameters are discontinuous, are all defined by two data lines (upper-side and
338lower-side values). This input format would also be needed for a graphic plot
339of the layered model. Layers which have different parameter values at top and
340bottom, will be treated as layers with a constant gradient, and will be
341discretised to a number of homogeneous sublayers. Errors due to the
342discretisation are limited within about 5%% (changeable, see peglobal.h).
343'''.lstrip()
345 return template % d
348class PoelError(Exception):
349 pass
352class PoelRunner(object):
354 def __init__(self, tmp=None):
356 self.tempdir = mkdtemp(prefix='poelrun', dir=tmp)
357 self.program = program_bins['poel']
358 self.config = None
360 def run(self, config):
361 self.config = config
363 input_fn = pjoin(self.tempdir, 'input')
365 f = open(input_fn, 'w')
366 poel_input = config.string_for_config()
368 logger.debug(
369 '===== begin poel input =====\n%s===== end poel input ====='
370 % poel_input)
372 f.write(poel_input)
373 f.close()
374 program = self.program
376 old_wd = os.getcwd()
378 os.chdir(self.tempdir)
380 try:
381 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
382 except OSError:
383 os.chdir(old_wd)
384 raise PoelError(
385 '''could not start poel executable: "%s"
386Available fomosto backends and download links to the modelling codes are listed
387on
389 https://pyrocko.org/docs/current/apps/fomosto/backends.html
391''' % program)
393 (poel_output, poel_error) = proc.communicate('input\n')
395 logger.debug(
396 '===== begin poel output =====\n%s===== end poel output ====='
397 % poel_output)
398 if poel_error:
399 logger.error(
400 '===== begin poel error =====\n%s===== end poel error ====='
401 % poel_error)
403 errmess = []
404 if proc.returncode != 0:
405 errmess.append(
406 'poel had a non-zero exit state: %i' % proc.returncode)
407 if poel_error:
408 errmess.append('poel emitted something via stderr')
409 if poel_output.lower().find('error') != -1:
410 errmess.append("the string 'error' appeared in poel output")
412 if errmess:
413 os.chdir(old_wd)
414 raise PoelError(
415 '''===== begin poel input =====\n%s===== end poel input =====
416===== begin poel output =====\n%s===== end poel output =====
417===== begin poel error =====\n%s===== end poel error =====
418%s
419poel has been invoked as "%s"''' % (
420 poel_input,
421 poel_output,
422 poel_error,
423 '\n'.join(errmess),
424 program))
426 self.poel_output = poel_output
427 self.poel_error = poel_error
429 os.chdir(old_wd)
431 def get_traces(self):
433 if self.config.sw_equidistant_x == 1:
434 nx = self.config.no_distances
435 xmin, xmax = self.config.distances
436 if nx > 1:
437 dx = (xmax-xmin)/(nx-1)
438 else:
439 dx = 1.0
440 distances = [xmin + ix*dx for ix in range(nx)]
441 else:
442 distances = self.config.distances
444 if self.config.sw_equidistant_z == 1:
445 nrz = self.config.no_depths
446 rzmin, rzmax = self.config.depths
447 if nrz > 1:
448 drz = (rzmax-rzmin)/(nrz-1)
449 else:
450 drz = 1.0
451 rdepths = [rzmin + irz*drz for irz in range(nrz)]
452 else:
453 rdepths = self.config.depths
455 sz = self.config.s_start_depth
457 fns = self.config.get_output_filenames(self.tempdir)
458 traces = []
459 for comp, fn in zip(poel_components, fns):
460 data = num.loadtxt(fn, skiprows=1, dtype=float)
461 nsamples, ntraces = data.shape
462 ntraces -= 1
463 tmin = data[0, 0]
464 deltat = (data[-1, 0] - data[0, 0])/(nsamples-1)
465 for itrace in range(ntraces):
466 x = distances[itrace % len(distances)]
467 rz = rdepths[itrace // len(distances)]
468 tr = trace.Trace(
469 '', '%i' % itrace, 'c', comp,
470 tmin=tmin, deltat=deltat, ydata=data[:, itrace+1],
471 meta={'itrace': itrace, 'x': x, 'rz': rz, 'sz': sz})
473 traces.append(tr)
475 return traces
477 def __del__(self):
478 shutil.rmtree(self.tempdir)
481class PoelGFBuilder(gf.builder.Builder):
482 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
483 force=False):
485 if block_size is None:
486 block_size = (51, 1, 51)
488 self.store = gf.store.Store(store_dir, 'w')
490 gf.builder.Builder.__init__(
491 self, self.store.config, step, block_size=block_size, force=force)
493 self.poel_config = self.store.get_extra('poel')
495 self.tmp = tmp
496 if self.tmp is not None:
497 util.ensuredir(self.tmp)
499 def cleanup(self):
500 self.store.close()
502 def work_block(self, index):
503 logger.info('Starting block %i / %i' % (index+1, self.nblocks))
505 runner = PoelRunner(tmp=self.tmp)
507 conf = PoelConfigFull(**self.poel_config.items())
509 (firstrz, sz, firstx), (lastrz, sz, lastx), (nrz, _, nx) = \
510 self.get_block_extents(index)
512 conf.s_start_depth = sz
513 conf.s_end_depth = sz
514 conf.sw_equidistant_x = 1
515 conf.distances = [firstx, lastx]
516 conf.sw_equidistant_z = 1
517 conf.no_distances = nx
518 conf.depths = [firstrz, lastrz]
519 conf.no_depths = nrz
520 conf.no_t_samples = int(
521 round(conf.t_window * self.gf_config.sample_rate)) + 1
523 runner.run(conf)
525 comp2ig = dict([(c, ig) for (ig, c) in enumerate(poel_components)])
527 rawtraces = runner.get_traces()
529 self.store.lock()
531 for tr in rawtraces:
533 x = tr.meta['x']
534 rz = tr.meta['rz']
535 sz = tr.meta['sz']
537 ig = comp2ig[tr.channel]
539 gf_tr = gf.store.GFTrace(
540 tr.get_ydata(),
541 int(round(tr.tmin / tr.deltat)),
542 tr.deltat)
544 self.store.put((rz, sz, x, ig), gf_tr)
546 self.store.unlock()
548 logger.info('Done with block %i / %i' % (index+1, self.nblocks))
551def init(store_dir, variant):
552 if variant is not None:
553 raise gf.store.StoreError('unsupported variant: %s' % variant)
555 poel = PoelConfig()
557 store_id = os.path.basename(os.path.realpath(store_dir))
559 config = gf.meta.ConfigTypeB(
560 modelling_code_id='poel',
561 id=store_id,
562 ncomponents=10,
563 component_scheme='poroelastic10',
564 sample_rate=0.1,
565 distance_min=10.0,
566 distance_max=20.0,
567 distance_delta=5.0,
568 source_depth_min=10.0,
569 source_depth_max=20.0,
570 source_depth_delta=5.0,
571 receiver_depth_min=10.0,
572 receiver_depth_max=10.0,
573 receiver_depth_delta=5.0)
575 return gf.store.Store.create_editables(
576 store_dir,
577 config=config,
578 extra={'poel': poel})
581def build(
582 store_dir,
583 force=False,
584 nworkers=None,
585 continue_=False,
586 step=None,
587 iblock=None):
589 return PoelGFBuilder.build(
590 store_dir,
591 force=force,
592 nworkers=nworkers,
593 continue_=continue_,
594 step=step,
595 iblock=iblock)