/[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 4576 - (show annotations)
Mon Dec 9 23:35:30 2013 UTC (5 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 3943 byte(s)
python3ified things, replaced mixed whitespace and xrange calls
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2013 by University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development since 2012 by School of Earth Sciences
11 #
12 ##############################################################################
13
14 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
15 http://www.uq.edu.au
16 Primary Business: Queensland, Australia"""
17 __license__="""Licensed under the Open Software License version 3.0
18 http://www.opensource.org/licenses/osl-3.0.php"""
19 __url__="https://launchpad.net/escript-finley"
20
21 from esys.escript import *
22 from esys.escript import unitsSI as U
23 from esys.escript.pdetools import Locator
24 from esys.finley import Brick, Rectangle
25 from esys.weipa import saveSilo
26 from esys.downunder import Ricker, SonicWave, SimpleSEGYWriter
27 from math import ceil
28
29
30 DIM=2 # spatial dimension
31
32 depth=1*U.km # depth
33 v_p_top=1.5*U.km/U.sec
34 v_p_bottom=3*U.km/U.sec
35 absorption_zone=300*U.m
36 ne_z=400
37
38 reflector_at=0.5*depth
39
40
41 t_end=1.*U.sec
42 frq=20.*U.Hz
43 sampling_interval=4*U.msec
44 numRcvPerLine=101
45 rangeRcv=800*U.m
46
47 # location of source in crossing array lines with in 0..numRcvInLine one needs to be None
48 srcEW=numRcvPerLine/2
49 srcNS=None
50
51 # dommain dimension
52 width_x=rangeRcv + 4*absorption_zone
53 width_y=width_x
54 #
55 # create array
56 #
57 receiver_line=[2*absorption_zone + i * (rangeRcv/(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
58 #
59 # set source location with tag "source""
60 #
61 src_tags=["source"]
62 if DIM == 2:
63 src_locations = [ (receiver_line[srcEW], depth)]
64 src_loc_2D=(receiver_line[srcEW], 0.)
65 else:
66 if srcEW:
67 srcNS=numRcvPerLine/2
68 elif srcNS:
69 srcNS=numRcvPerLine/2
70 else:
71 raise ValueError("on of the variables srcEW or srcNS must be None!")
72 src_locations = [ (receiver_line[srcEW], receiver_line[srcNS], depth)]
73 src_loc_2D=(receiver_line[srcEW], receiver_line[srcNS])
74 #
75 # create sensor arrays:
76 #
77 # East-west line of receiver
78 rcvEW_locations=[]
79 rgEW=[]
80 mid_point=receiver_line[len(receiver_line)/2]
81
82 for ix in range(len(receiver_line)):
83 if DIM == 2:
84 rcvEW_locations.append((receiver_line[ix], depth))
85 rgEW.append( ( receiver_line[ix], 0.) )
86 else:
87 rcvEW_locations.append((receiver_line[ix], mid_point, depth))
88 rgEW.append( ( receiver_line[ix], mid_point) )
89 # North-south line of receiver
90 if DIM == 3:
91 rcvNS_locations=[]
92 rgNS=[]
93
94 for iy in range(len(receiver_line)):
95 rcvNS_locations.append((mid_point, receiver_line[iy], depth))
96 rgNS.append( ( mid_point, receiver_line[iy]) )
97 #
98 # create domain:
99 #
100 if DIM == 2:
101 domain=Rectangle(ceil(ne_z*width_x/depth),ne_z,l0=width_x,l1=depth,
102 diracPoints=src_locations, diracTags=src_tags)
103 else:
104 domain=Brick(ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),ne_z,l0=width_x,l1=width_y,l2=depth,
105 diracPoints=src_locations, diracTags=src_tags)
106 wl=Ricker(frq)
107 m=whereNegative(Function(domain).getX()[DIM-1]-reflector_at)
108 v_p=v_p_bottom*m+v_p_top*(1-m)
109
110 sw=SonicWave(domain, v_p, source_tag=src_tags[0], wavelet=wl, absorption_zone=absorption_zone, lumping=True)
111
112 locEW=Locator(domain,rcvEW_locations)
113 tracerEW=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval)
114 if DIM==3:
115 locNS=Locator(domain,rcvNS_locations)
116 tracerNS=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D, sampling_interval=sampling_interval)
117
118
119 t=0.
120 mkDir('tmp')
121 n=0
122 while t < t_end:
123 t,p = sw.update(t+sampling_interval)
124 tracerEW.addRecord(locEW(p))
125 if DIM==3: tracerNS.addRecord(locNS(p))
126 print t, locEW(p)[:4], wl.getValue(t)
127 if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n/5,), p=p)
128 n+=1
129 tracerEW.write('lineEW.sgy')
130 if DIM == 3: tracerNS.write('lineNS.sgy')

  ViewVC Help
Powered by ViewVC 1.1.26