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

Contents of /trunk/doc/examples/inversion/synthetic_sonicHTI.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: 6371 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 from __future__ import division
2 from __future__ import print_function
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2018 by The University of Queensland
6 # http://www.uq.edu.au
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Apache License, version 2.0
10 # http://www.apache.org/licenses/LICENSE-2.0
11 #
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16 from __future__ import print_function
17
18 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Apache License, version 2.0
22 http://www.apache.org/licenses/LICENSE-2.0"""
23 __url__="https://launchpad.net/escript-finley"
24
25 from esys.escript import *
26 from esys.escript import unitsSI as U
27 from esys.escript.pdetools import Locator
28 from esys.weipa import saveSilo
29 from esys.downunder import Ricker, SonicHTIWave, SimpleSEGYWriter
30 from math import ceil
31 import time, os
32
33 try:
34 from esys.ripley import Brick, Rectangle
35 HAVE_RIPLEY = True
36 except ImportError:
37 HAVE_RIPLEY = False
38
39 if HAVE_RIPLEY:
40 DIM=2 # spatial dimension
41
42 # layers from the bottom up:
43 layers = [ 1*U.km , 1*U.km ,700*U.m, 500*U.m, 800*U.m ]
44 v_Ps= [ 3.8 * U.km/U.sec , 3. * U.km/U.sec, 2.5*U.km/U.sec, 1.9*U.km/U.sec, 1.5*U.km/U.sec]
45 epss =[ 0., 0.24, 0, 0.1, 0]
46 deltas=[ 0., 0.1, 0.,0.03,0 ]
47 azmths=[ 0.,0.,0, 0, 0.]
48
49 dt=0.5*U.msec
50
51 ne_z=40
52
53 dt=0.5*U.msec
54
55 t_end=0.008*U.sec #only this low for testing purposes
56 frq=15.*U.Hz
57 tcenter=None
58 sampling_interval=4*U.msec
59 numRcvPerLine=101
60 rangeRcv=4.*U.km
61 src_dir=[0,1]
62 absorption_zone=1000*U.m
63
64 # location of source in crossing array lines with in 0..numRcvInLine one needs to be None
65 srcEW=numRcvPerLine//2
66 srcNS=None
67 # dommain dimension
68 width_x=rangeRcv + 2*absorption_zone
69 width_y=width_x
70 depth=sum(layers)
71 ne_x=int(ceil(ne_z*width_x/depth))
72 #
73 # create array
74 #
75 receiver_line=[ absorption_zone + i * (rangeRcv//(numRcvPerLine-1) ) for i in range(numRcvPerLine) ]
76 #
77 # set source location with tag "source""
78 #
79 src_tags=["source"]
80
81 if srcEW:
82 srcNS=numRcvPerLine//2
83 elif srcNS:
84 srcEW=numRcvPerLine//2
85 else:
86 raise ValueError("on of the variables srcEW or srcNS must be None!")
87 if DIM == 2:
88 src_locations = [ (receiver_line[srcEW], depth) ]
89 src_loc_2D=(receiver_line[srcEW], 0.)
90 else:
91 src_locations = [ (receiver_line[srcEW], receiver_line[srcNS], depth)]
92 src_loc_2D=(receiver_line[srcEW], receiver_line[srcNS])
93
94 #
95 # create sensor arrays:
96 #
97 # East-west line of receiver
98 rcv_locations=[]
99 rg=[]
100 mid_point=receiver_line[len(receiver_line)//2]
101
102 for ix in range(len(receiver_line)):
103 if DIM == 2:
104 rcv_locations.append((receiver_line[ix], depth))
105 rg.append( ( receiver_line[ix], 0.) )
106 else:
107 rcv_locations.append((receiver_line[ix], mid_point, depth))
108 rg.append( ( receiver_line[ix], mid_point) )
109 # North-south line of receiver
110 if DIM == 3:
111 for iy in range(len(receiver_line)):
112 rcv_locations.append((mid_point, receiver_line[iy], depth))
113 rg.append( ( mid_point, receiver_line[iy]) )
114 #
115 # create domain:
116 #
117 if DIM == 2:
118 domain=Rectangle(ne_x, ne_z ,l0=width_x, l1=depth,
119 diracPoints=src_locations, diracTags=src_tags)
120 else:
121 domain=Brick(ne_x,ne_x,ne_z,l0=width_x,l1=width_y,l2=depth,
122 diracPoints=src_locations, diracTags=src_tags)
123 wl=Ricker(frq, tcenter)
124
125 #======================================================================
126 z=Function(domain).getX()[DIM-1]
127 z_bottom=0
128 v_p=0
129 delta=0
130 vareps=0
131 azmth=0
132 rho=0
133 for l in range(len(layers)):
134 m=wherePositive(z-z_bottom)*whereNonPositive(z-(z_bottom+layers[l]))
135 v_p=v_p*(1-m)+v_Ps[l]*m
136 vareps=vareps*(1-m)+epss[l]*m
137 azmth=azmth*(1-m)+azmths[l]*m
138 delta=delta*(1-m)+deltas[l]*m
139 z_bottom+=layers[l]
140
141 sw=SonicHTIWave(domain, v_p, wl, src_tags[0], dt=dt, source_vector = src_dir, eps=vareps, delta=delta, azimuth=azmth, \
142 absorption_zone=absorption_zone, absorption_cut=1e-2, lumping=False)
143
144 #
145 # print some info:
146 #
147 print("ne_x = ", ne_x)
148 print("ne_z = ", ne_z)
149 print("h_x = ", width_x/ne_x)
150 print("h_z = ", depth/ne_z)
151 print("dt = ", sw.getTimeStepSize()*1000, "msec")
152 print("width_x = ", width_x)
153 print("depth = ", depth)
154 print("number receivers = ", numRcvPerLine)
155 print("receiver spacing = ", receiver_line[1]-receiver_line[0])
156 print("sampling time = ", sampling_interval*1000,"msec")
157 print("source @ ", src_locations[0])
158 #
159 loc=Locator(domain,rcv_locations)
160 tracerP=SimpleSEGYWriter(receiver_group=rg, source=src_loc_2D, sampling_interval=sampling_interval, text='P')
161 tracerQ=SimpleSEGYWriter(receiver_group=rg, source=src_loc_2D, sampling_interval=sampling_interval, text='Q')
162
163 if not tracerP.obspy_available():
164 print("\nWARNING: obspy not available, SEGY files will not be written\n")
165 elif getMPISizeWorld() > 1:
166 print("\nWARNING: SEGY files cannot be written with multiple processes\n")
167
168 t=0.
169 OUT_DIR="out%sm%smus"%(int(width_x/ne_x),int(sw.getTimeStepSize()*1000000))
170 mkDir(OUT_DIR)
171 n=0
172 k=0
173 timer1=time.time()
174 while t < t_end:
175 t,u = sw.update(t+sampling_interval)
176 Plog=loc(u[1])
177 Qlog=loc(u[0])
178 tracerP.addRecord(Plog)
179 tracerQ.addRecord(Qlog)
180 print(t, wl.getValue(t)," :", Plog[0], Plog[srcEW], Plog[-1])
181 timer1=time.time()-timer1
182 print("time= %e sec; %s sec per step"%(timer1,timer1/max(sw.n,1)))
183
184 if tracerP.obspy_available() and getMPISizeWorld() == 1:
185 tracerP.write(os.path.join(OUT_DIR,'lineP.sgy'))
186 tracerQ.write(os.path.join(OUT_DIR,'lineQ.sgy'))
187
188 else: # no ripley
189 print("The Ripley module is not available")
190

  ViewVC Help
Powered by ViewVC 1.1.26