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

Contents of /trunk/doc/examples/inversion/synthetic_VTI.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: 6677 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-2018 by The 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, VTIWave, SimpleSEGYWriter
28 from math import ceil
29
30 try:
31 from esys.speckley import Brick, Rectangle
32 HAVE_SPECKLEY=True
33 except ImportError:
34 HAVE_SPECKLEY=False
35
36 if HAVE_SPECKLEY:
37 DIM=2 # spatial dimension
38
39 depth=1*U.km # depth
40 v_p_top=1.5*U.km/U.sec
41 v_p_bottom=3*U.km/U.sec
42 absorption_zone=100*U.m
43 ne_z=50.
44
45 reflector_at=0.5*depth
46
47
48 t_end=0.008*U.sec #only this low for testing purposes
49 frq=8.*U.Hz
50 sampling_interval=4*U.msec
51 numRcvPerLine=101
52 rangeRcv=800*U.m
53
54 # location of source in crossing array lines with in 0..numRcvInLine one needs to be None
55 srcEW=numRcvPerLine//2
56 srcNS=None
57
58 # dommain dimension
59 width_x=rangeRcv + 4*absorption_zone
60 width_y=width_x
61 #
62 # create array
63 #
64 receiver_line=[2*absorption_zone + i * (rangeRcv//(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
65 #
66 # set source location with tag "source""
67 #
68 src_tags=["source"]
69 if DIM == 2:
70 src_locations = [ (receiver_line[srcEW], depth)]
71 src_loc_2D=(receiver_line[srcEW], 0.)
72 else:
73 if srcEW:
74 srcNS=numRcvPerLine//2
75 elif srcNS:
76 srcEW=numRcvPerLine//2
77 else:
78 raise ValueError("on of the variables srcEW or srcNS must be None!")
79 src_locations = [ (receiver_line[srcEW], receiver_line[srcNS], depth)]
80 src_loc_2D=(receiver_line[srcEW], receiver_line[srcNS])
81 #
82 # create sensor arrays:
83 #
84 # East-west line of receiver
85 rcvEW_locations=[]
86 rgEW=[]
87 mid_point=receiver_line[len(receiver_line)//2]
88
89 for ix in range(len(receiver_line)):
90 if DIM == 2:
91 rcvEW_locations.append((receiver_line[ix], depth))
92 rgEW.append( ( receiver_line[ix], 0.) )
93 else:
94 rcvEW_locations.append((receiver_line[ix], mid_point, depth))
95 rgEW.append( ( receiver_line[ix], mid_point) )
96 # North-south line of receiver
97 if DIM == 3:
98 rcvNS_locations=[]
99 rgNS=[]
100
101 for iy in range(len(receiver_line)):
102 rcvNS_locations.append((mid_point, receiver_line[iy], depth))
103 rgNS.append( ( mid_point, receiver_line[iy]) )
104 #
105 # create domain:
106 #
107 order = 5
108 if DIM == 2:
109 domain=Rectangle(order, ceil(ne_z*width_x/depth),ne_z,l0=width_x,l1=depth,
110 diracPoints=src_locations, diracTags=src_tags)
111 else:
112 domain=Brick(order, ceil(ne_z*width_x/depth), ceil(ne_z*width_y/depth),
113 ne_z, l0=width_x, l1=width_y, l2=depth,
114 diracPoints=src_locations, diracTags=src_tags)
115 wl=Ricker(frq)
116
117 #======================================================================
118 # m=whereNegative(Function(domain).getX()[DIM-1]-reflector_at)
119 # v_p=v_p_bottom*m+v_p_top*(1-m)
120 v_p=2*U.km/U.sec
121 v_s=0.9*U.km/U.sec
122 vareps=0.1*0
123 gamma=0.15*0
124 delta=0.05*0
125 rho=2000*U.kg/U.m**3
126 src_dir=[0,0,1]
127
128 sw=VTIWave(domain, v_p, v_s, wl, src_tags[0], source_vector = src_dir,
129 eps=vareps, gamma=gamma, delta=delta, rho=rho,
130 absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True)
131
132 locEW=Locator(domain,rcvEW_locations)
133 tracerEW_x=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='x-displacement - east-west line')
134 tracerEW_z=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='z-displacement - east-west line')
135 if DIM==3:
136 locNS=Locator(domain,rcvNS_locations)
137 tracerEW_y=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D,
138 sampling_interval=sampling_interval,
139 text='x-displacement - east-west line')
140 tracerNS_x=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
141 sampling_interval=sampling_interval,
142 text='x-displacement - north-south line')
143 tracerNS_y=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
144 sampling_interval=sampling_interval,
145 text='y-displacement - north-south line')
146 tracerNS_z=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D,
147 sampling_interval=sampling_interval,
148 text='z-displacement - north-south line')
149 if not tracerEW_x.obspy_available():
150 print("\nWARNING: obspy not available, SEGY files will not be written\n")
151 elif getMPISizeWorld() > 1:
152 print("\nWARNING: SEGY files cannot be written with multiple processes\n")
153
154 t=0.
155 mkDir('tmp')
156 n=0
157 while t < t_end:
158 t,u = sw.update(t+sampling_interval)
159 tracerEW_x.addRecord(locEW(u[0]))
160 tracerEW_z.addRecord(locEW(u[DIM-1]))
161 if DIM==3:
162 tracerEW_y.addRecord(locEW(u[1]))
163 tracerNS_x.addRecord(locNS(u[0]))
164 tracerNS_y.addRecord(locNS(u[1]))
165 tracerNS_z.addRecord(locNS(u[2]))
166 print(t, locEW(u[DIM-1])[len(rgEW)//2-4:len(rgEW)//2+1], wl.getValue(t))
167 #if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n/5,), u=u)
168 saveSilo("tmp/u_%d.silo"%(n,), u=u, cycle=n, time=t)
169 n+=1
170 if tracerEW_x.obspy_available() and getMPISizeWorld() == 1:
171 tracerEW_x.write('lineEW_x.sgy')
172 tracerEW_z.write('lineEW_z.sgy')
173 if DIM == 3:
174 tracerEW_y.write('lineEW_y.sgy')
175 tracerNS_x.write('lineNS_x.sgy')
176 tracerNS_y.write('lineNS_y.sgy')
177 tracerNS_z.write('lineNS_z.sgy')
178
179 else: # no speckley
180 print("The Speckley module is not available")
181

  ViewVC Help
Powered by ViewVC 1.1.26