/[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 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: 7546 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-2015 by 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, SimpleSEGYWriter, HTIWave
28 from math import ceil
29 from time import time
30
31 try:
32 from esys.speckley import Brick, Rectangle
33 HAVE_SPECKLEY=True
34 except ImportError:
35 HAVE_SPECKLEY=False
36
37 if HAVE_SPECKLEY:
38 DIM=2 # spatial dimension
39
40
41 ORDER = 5
42 ne_z= 20
43
44 # layers from the bottom up:
45 layers=[20*U.m, 180*U.m ]
46 v_Ps=[i*U.km/U.sec for i in [3, 2.5]]
47 v_Ss= [i*U.km/U.sec for i in [3, 2]]
48 rhos=[i*U.kg/U.m**3 for i in [2.6, 2.1]]
49 epss=[0, .110]
50 gammas=[0, 0.035]
51 deltas=[0, 0.255]
52 src_dir=[0,0,1]
53
54 t_end=0.01*U.sec #increase this end time as desired
55 frq=50.*U.Hz
56 sampling_interval=2*U.msec
57 numRcvPerLine=101
58 rangeRcv=200*U.m
59
60
61 # location of source
62 if DIM == 2:
63 src_locations = [(0, 0)]
64 else:
65 src_locations = [(0, 0, 0)]
66
67 # domain dimensions
68 width_x=rangeRcv
69 width_y=width_x
70 depth=sum(layers)
71 #
72 # create array
73 #
74 receiver_line=[i * (rangeRcv/(numRcvPerLine-1)) for i in range(numRcvPerLine)]
75 #
76 # set source location with tag "source""
77 #
78 src_tags=["source"]
79
80 src_loc_2D=(0, 0)
81
82
83
84 #
85 # create sensor arrays:
86 #
87 # East-west line of receivers
88 rcvEW_locations=[]
89 # North-south line of receivers (if 3 dimensional problem)
90 rcvNS_locations=[]
91 rgEW=[]
92 rgNS=[]
93 mid_point=receiver_line[len(receiver_line)//2]
94
95 for ix in range(len(receiver_line)):
96 rgEW.append((receiver_line[ix], 0))
97 if DIM == 2:
98 rcvEW_locations.append((receiver_line[ix], 0))
99 else:
100 rcvEW_locations.append((receiver_line[ix], 0, 0))
101 rcvNS_locations.append((0, receiver_line[ix], 0))
102 rgNS.append((0, receiver_line[ix]))
103 # North-south line of receivers
104 if DIM == 3:
105 for iy in range(len(receiver_line)):
106 rcv_locations.append((mid_point, receiver_line[iy], 0))
107 rg.append( ( mid_point, receiver_line[iy]) )
108 #
109 # create domain:
110 #
111 if DIM == 2:
112 domain = Rectangle(ORDER,
113 ceil(ne_z*width_x/depth), ne_z ,l0=width_x, l1=(-depth,0),
114 diracPoints=src_locations, diracTags=src_tags)
115 #suppress the x-component on the x boundary
116 q = whereZero(domain.getX()[0])*[1,0]
117 else:
118 domain=Brick(ORDER,
119 ceil(ne_z*width_x/depth), ceil(ne_z*width_y/depth), ne_z,
120 l0=width_x, l1=width_y, l2=(-depth,0),
121 diracPoints=src_locations, diracTags=src_tags)
122 q = wherePositive(
123 #suppress the x-component on the x boundary
124 whereZero(domain.getX()[0])*[1,0,0]
125 + #logical or
126 #suppress the y-component on the y boundary at the source
127 whereZero(domain.getX()[1])*[0,1,0])
128
129 # set up reciever locations
130 locEW=Locator(domain,rcvEW_locations)
131 tracerEW_x=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D,
132 sampling_interval=sampling_interval,
133 text='x-displacement - east-west line')
134 tracerEW_z=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D,
135 sampling_interval=sampling_interval,
136 text='z-displacement - east-west line')
137 if DIM==3:
138 locNS=Locator(domain,rcvNS_locations)
139 tracerEW_y=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D,
140 sampling_interval=sampling_interval,
141 text='x-displacement - east-west line')
142 tracerNS_x=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
143 sampling_interval=sampling_interval,
144 text='x-displacement - north-south line')
145 tracerNS_y=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
146 sampling_interval=sampling_interval,
147 text='y-displacement - north-south line')
148 tracerNS_z=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
149 sampling_interval=sampling_interval,
150 text='z-displacement - north-south line')
151 if not tracerEW_x.obspy_available():
152 print("\nWARNING: obspy not available, SEGY files will not be written\n")
153 elif getMPISizeWorld() > 1:
154 print("\nWARNING: SEGY files cannot be written with multiple processes\n")
155
156
157 #======================================================================
158 z=ReducedFunction(domain).getX()[DIM-1]
159 z_bottom=-depth
160 v_p=0
161 v_s=0
162 delta=0
163 vareps=0
164 gamma=0
165 rho=0
166 for l in range(len(layers)):
167 m=wherePositive(z-z_bottom)*whereNonPositive(z-(z_bottom+layers[l]))
168 v_p=v_p*(1-m)+v_Ps[l]*m
169 v_s=v_s*(1-m)+v_Ss[l]*m
170 rho=rho*(1-m)+rhos[l]*m
171 vareps=vareps*(1-m)+epss[l]*m
172 gamma=gamma*(1-m)+gammas[l]*m
173 delta=delta*(1-m)+deltas[l]*m
174 z_bottom+=layers[l]
175
176 wl=Ricker(frq)
177 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wl.getTimeScale())
178
179 sw=HTIWave(domain, v_p, v_s, wl, src_tags[0], source_vector = src_dir,
180 eps=vareps, gamma=gamma, delta=delta, rho=rho,
181 absorption_zone=None, absorption_cut=1e-2, lumping=True, dt=dt)
182 sw.setQ(q)
183
184 locEW=Locator(domain, rcvEW_locations)
185 if DIM == 3:
186 locNS=Locator(domain, rcvNS_locations)
187
188 mkDir('output')
189
190 t=0.
191 n=0
192 k=0
193 u=None
194 while t < t_end:
195 start = time()
196 t,u = sw.update(t+sampling_interval)
197 tracerEW_x.addRecord(locEW(u[0]))
198 tracerEW_z.addRecord(locEW(u[DIM-1]))
199 if DIM==3:
200 tracerEW_y.addRecord(locEW(u[1]))
201 tracerNS_x.addRecord(locNS(u[0]))
202 tracerNS_y.addRecord(locNS(u[1]))
203 tracerNS_z.addRecord(locNS(u[2]))
204 print(t, locEW(u[DIM-1])[len(rgEW)//2-4:len(rgEW)//2+1], wl.getValue(t))
205 k+=1
206 if k%5 == 0:
207 saveSilo("output/normalHTI_%d.silo"%(n,), v_p=v_p, u=u, cycle=k, time=t)
208 n += 1
209 if k%5 != 0:
210 saveSilo("output/normalHTI_%d.silo"%(n,), v_p=v_p, u=u, cycle=k, time=t)
211 if tracerEW_x.obspy_available() and getMPISizeWorld() == 1:
212 tracerEW_x.write('output/lineEW_x.sgy')
213 tracerEW_z.write('output/lineEW_z.sgy')
214 if DIM == 3:
215 tracerEW_y.write('output/lineEW_y.sgy')
216 tracerNS_x.write('output/lineNS_x.sgy')
217 tracerNS_y.write('output/lineNS_y.sgy')
218 tracerNS_z.write('output/lineNS_z.sgy')
219
220 else: # no speckley
221 print("The Speckley module is not available")
222

  ViewVC Help
Powered by ViewVC 1.1.26