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