/[escript]/trunk/doc/examples/inversion/synthetic_sonic.py
ViewVC logotype

Contents of /trunk/doc/examples/inversion/synthetic_sonic.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 4961 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 ##############################################################################
2 #
3 # Copyright (c) 2003-2018 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Apache License, version 2.0
8 # http://www.apache.org/licenses/LICENSE-2.0
9 #
10 # Development 2012-2013 by School of Earth Sciences
11 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
12 #
13 ##############################################################################
14 from __future__ import print_function, division
15
16 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Apache License, version 2.0
20 http://www.apache.org/licenses/LICENSE-2.0"""
21 __url__="https://launchpad.net/escript-finley"
22
23 from esys.escript import *
24 from esys.escript import unitsSI as U
25 from esys.escript.pdetools import Locator
26 from esys.weipa import saveSilo
27 from esys.downunder import Ricker, SonicWave, SimpleSEGYWriter
28 from math import ceil
29 try:
30 from esys.speckley import Brick, Rectangle
31 HAVE_SPECKLEY=True
32 except ImportError:
33 HAVE_SPECKLEY=False
34
35 if HAVE_SPECKLEY:
36 DIM=2 # spatial dimension
37
38 depth=1*U.km # depth
39 v_p_top=1.5*U.km/U.sec
40 v_p_bottom=3*U.km/U.sec
41 absorption_zone=300*U.m
42 ne_z=40
43
44 reflector_at=0.5*depth
45
46
47 t_end=0.008*U.sec #only this low for testing purposes
48 frq=20.*U.Hz
49 sampling_interval=4*U.msec
50 numRcvPerLine=101
51 rangeRcv=800*U.m
52
53 # location of source in crossing array lines with in 0..numRcvInLine one needs to be None
54 srcEW=numRcvPerLine//2
55 srcNS=None
56
57 # dommain dimension
58 width_x=rangeRcv + 4*absorption_zone
59 width_y=width_x
60 #
61 # create array
62 #
63 receiver_line=[2*absorption_zone + i * (rangeRcv//(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
64 #
65 # set source location with tag "source""
66 #
67 src_tags=["source"]
68 if DIM == 2:
69 src_locations = [ (receiver_line[srcEW], depth)]
70 src_loc_2D=(receiver_line[srcEW], 0.)
71 else:
72 if srcEW:
73 srcNS=numRcvPerLine//2
74 elif srcNS:
75 srcNS=numRcvPerLine//2
76 else:
77 raise ValueError("on of the variables srcEW or srcNS must be None!")
78 src_locations = [ (receiver_line[srcEW], receiver_line[srcNS], depth)]
79 src_loc_2D=(receiver_line[srcEW], receiver_line[srcNS])
80 #
81 # create sensor arrays:
82 #
83 # East-west line of receiver
84 rcvEW_locations=[]
85 rgEW=[]
86 mid_point=receiver_line[len(receiver_line)//2]
87
88 for ix in range(len(receiver_line)):
89 if DIM == 2:
90 rcvEW_locations.append((receiver_line[ix], depth))
91 rgEW.append( ( receiver_line[ix], 0.) )
92 else:
93 rcvEW_locations.append((receiver_line[ix], mid_point, depth))
94 rgEW.append( ( receiver_line[ix], mid_point) )
95 # North-south line of receiver
96 if DIM == 3:
97 rcvNS_locations=[]
98 rgNS=[]
99
100 for iy in range(len(receiver_line)):
101 rcvNS_locations.append((mid_point, receiver_line[iy], depth))
102 rgNS.append( ( mid_point, receiver_line[iy]) )
103 #
104 # create domain:
105 #
106 order = 5
107 if DIM == 2:
108 domain=Rectangle(order, ceil(ne_z*width_x/depth),ne_z,l0=width_x,l1=depth,
109 diracPoints=src_locations, diracTags=src_tags)
110 else:
111 domain=Brick(order, ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),ne_z,l0=width_x,l1=width_y,l2=depth,
112 diracPoints=src_locations, diracTags=src_tags)
113 wl=Ricker(frq)
114 m=whereNegative(Function(domain).getX()[DIM-1]-reflector_at)
115 v_p=v_p_bottom*m+v_p_top*(1-m)
116
117 sw=SonicWave(domain, v_p, source_tag=src_tags[0], wavelet=wl, absorption_zone=absorption_zone, lumping=True)
118
119 locEW=Locator(domain,rcvEW_locations)
120 tracerEW=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval)
121 if DIM==3:
122 locNS=Locator(domain,rcvNS_locations)
123 tracerNS=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D, sampling_interval=sampling_interval)
124
125 if not tracerEW.obspy_available():
126 print("\nWARNING: obspy not available, SEGY files will not be written\n")
127 elif getMPISizeWorld() > 1:
128 print("\nWARNING: SEGY files cannot be written with multiple processes\n")
129
130 t=0.
131 mkDir('tmp')
132 n=0
133 while t < t_end:
134 t,p = sw.update(t+sampling_interval)
135 tracerEW.addRecord(locEW(p))
136 if DIM==3: tracerNS.addRecord(locNS(p))
137 print(t, locEW(p)[:4], wl.getValue(t))
138 if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n//5,), p=p)
139 n+=1
140 if tracerEW.obspy_available() and getMPISizeWorld() == 1:
141 tracerEW.write('lineEW.sgy')
142 if DIM == 3:
143 tracerNS.write('lineNS.sgy')
144
145 else: # no speckley
146 print("The Speckley module is not available")
147

  ViewVC Help
Powered by ViewVC 1.1.26