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

Contents of /trunk/doc/examples/inversion/synthetic_HTI.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: 5433 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, HTIWave, SimpleSEGYWriter
27 from math import ceil
28
29
30 DIM=3 # spatial dimension
31
32 v_p_top=1.5*U.km/U.sec
33 v_p_bottom=3*U.km/U.sec
34 absorption_zone=300*U.m
35 ne_z=500.
36
37 if DIM==3:
38 # layers from the bottom up:
39 layers=[500*U.m, 500 * U.m ]
40 v_Ps=[ 3*U.km/U.sec, 2.*U.km/U.sec]
41 v_Ss=[1.3*U.km/U.sec, 0.9*U.km/U.sec ]
42 rhos=[2000*U.kg/U.m**3, 2000*U.kg/U.m**3]
43 epss=[0., 0.1]
44 gammas=[0., 0.03]
45 deltas=[0., 0.1]
46
47 depth=sum(layers)
48
49 src_dir=[0,0,1]
50 else:
51 # this is for DIM=2 case:
52 v_p=2*U.km/U.sec
53 v_s=0.9*U.km/U.sec
54 rho=2000*U.kg/U.m**3
55
56 vareps=0.1
57 gamma=0.15
58 #delta=-0.1
59 delta=0.0
60 delta=0.1
61 depth=1*U.km
62 src_dir=[1,1,0]
63
64 t_end=0.8*U.sec
65 frq=20.*U.Hz
66 sampling_interval=4*U.msec
67 numRcvPerLine=101
68 rangeRcv=800*U.m
69
70
71 # location of source in crossing array lines with in 0..numRcvInLine one needs to be None
72 srcEW=numRcvPerLine/2
73 srcNS=None
74
75 # dommain dimension
76 width_x=rangeRcv + 4*absorption_zone
77 width_y=width_x
78
79 #
80 # create array
81 #
82 receiver_line=[2*absorption_zone + i * (rangeRcv/(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
83 #
84 # set source location with tag "source""
85 #
86 src_tags=["source"]
87
88 if srcEW:
89 srcNS=numRcvPerLine/2
90 elif srcNS:
91 srcEW=numRcvPerLine/2
92 else:
93 raise ValueError("on of the variables srcEW or srcNS must be None!")
94 src_loc_2D=(receiver_line[srcEW], receiver_line[srcNS])
95 if DIM == 2:
96 src_locations = [ src_loc_2D ]
97 else:
98 src_locations = [ (receiver_line[srcEW], receiver_line[srcNS], depth)]
99 #
100 # create sensor arrays:
101 #
102 # East-west line of receiver
103 rcv_locations=[]
104 rg=[]
105 mid_point=receiver_line[len(receiver_line)/2]
106
107 for ix in range(len(receiver_line)):
108 if DIM == 2:
109 rcv_locations.append((receiver_line[ix], mid_point))
110 rg.append( ( receiver_line[ix], mid_point) )
111 else:
112 rcv_locations.append((receiver_line[ix], mid_point, depth))
113 rg.append( ( receiver_line[ix], mid_point) )
114 # North-south line of receiver
115 for iy in range(len(receiver_line)):
116 if DIM == 2:
117 rcv_locations.append((mid_point, receiver_line[iy]))
118 rg.append( ( mid_point, receiver_line[iy]) )
119 else:
120 rcv_locations.append((mid_point, receiver_line[iy], depth))
121 rg.append( ( mid_point, receiver_line[iy]) )
122 #
123 # create domain:
124 #
125 if DIM == 2:
126 domain=Rectangle(ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),l0=width_x,l1=width_y,
127 diracPoints=src_locations, diracTags=src_tags)
128 else:
129 domain=Brick(ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),ne_z,l0=width_x,l1=width_y,l2=depth,
130 diracPoints=src_locations, diracTags=src_tags)
131 wl=Ricker(frq)
132
133 #======================================================================
134 if DIM == 3:
135 z=Function(domain).getX()[DIM-1]
136 z_bottom=0
137 v_p=0
138 delta=0
139 vareps=0
140 gamma=0
141 rho=0
142 for l in range(len(layers)):
143 m=wherePositive(z-z_bottom)*whereNonNegative(z-(z_bottom+layers[l]))
144 v_p=v_p*(1-m)+v_Ps[l]*m
145 v_s=v_s*(1-m)+v_Ss[l]*m
146 rho=rho*(1-m)+rhos[l]*m
147 vareps=vareps*(1-m)+epss[l]*m
148 gamma=gamma*(1-m)+gammas[l]*m
149 delta=delta*(1-m)+deltas[l]*m
150 z_bottom+=layers[l]
151
152 sw=HTIWave(domain, v_p, v_s, wl, src_tags[0], source_vector = src_dir, eps=vareps, gamma=gamma, delta=delta, rho=rho, \
153 absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True)
154
155 loc=Locator(domain,rcv_locations)
156 tracer_x=SimpleSEGYWriter(receiver_group=rg, source=src_loc_2D, sampling_interval=sampling_interval, text='x-displacement')
157 tracer_y=SimpleSEGYWriter(receiver_group=rg, source=src_loc_2D, sampling_interval=sampling_interval, text='y-displacement')
158 if DIM==3:
159 tracer_z=SimpleSEGYWriter(receiver_group=rg, source=src_loc_2D, sampling_interval=sampling_interval, text='z-displacement')
160
161 t=0.
162 mkDir('tmp')
163 n=0
164 k=0
165 while t < t_end:
166 t,u = sw.update(t+sampling_interval)
167 tracer_x.addRecord(loc(u[0]))
168 tracer_y.addRecord(loc(u[1]))
169 if DIM==3:
170 tracer_z.addRecord(loc(u[2]))
171 print t, loc(u[0])[len(rg)/2-4:len(rg)/2+1], wl.getValue(t)
172 #if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n/5,), u=u)
173 if t>0.3 and t< 0.5:
174 saveSilo("tmp/u_%d.silo"%(k,), u=u)
175 k+=1
176 n+=1
177 tracer_x.write('line_x.sgy')
178 tracer_y.write('line_y.sgy')
179 if DIM == 3:
180 tracer_z.write('line_z.sgy')
181

  ViewVC Help
Powered by ViewVC 1.1.26