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

Contents of /trunk/doc/examples/inversion/synthetic_TTI.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: 7020 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, TTIWave, SimpleSEGYWriter
28 from math import ceil
29 import time
30 try:
31 from esys.speckley import Rectangle
32 HAVE_SPECKLEY=True
33 except ImportError:
34 HAVE_SPECKLEY=False
35
36 if HAVE_SPECKLEY:
37 # these are the layers from the top down
38 layers = [ 400*U.m , 100*U.m , 1.*U.km, ]
39 v_P= [ 2.86* U.km/U.sec , 1.5 * U.km/U.sec, 2.86 * U.km/U.sec ]
40 v_S= [ 1.79 * U.km/U.sec , 0.7* U.km/U.sec, 1.8*U.km/U.sec ]
41 eps = [ 0. , 0.5, 0.1 ]
42 delta= [ 0. , 0.5 , 0. ]
43 tilt= [ 0. , 0. , 0. ]
44 rho= [ 2000 * U.kg/U.m**3 , 2000 * U.kg/U.m**3, 2000 * U.kg/U.m**3 ]
45 #
46 # other input:
47 #
48 t_end=0.008*U.sec # only this low for testing purposes
49 frq=10.*U.Hz # dominant frequnce in the Ricker (maximum frequence ~ 2 * frq)
50 sampling_interval=4*U.msec # sampling interval
51 ne_z=None # number of elements in vertical direction, if none it is guessed
52 n_out = 5 # a silo file is created every n_out's sample
53 absorption_zone=100*U.m # absorbtion zone to be added in horizontal direction to the area covered by receiver line
54 # and subtracted from the lowest layer.
55 # defines the receiver line
56 rangeRcv=800*U.m # width of the receiver line
57 numRcvPerLine=101 # total number of receivers
58 src_id=numRcvPerLine//2 # location of source in crossing array lines with in 0..numRcvInLine
59 lumping = True
60 src_dir=[0,1]
61
62 # domain dimension
63 width_x=rangeRcv + 4*absorption_zone
64 depth=sum(layers)
65 if ne_z is None:
66 ne_z=int(ceil(depth*(2*frq)/min(v_P)))
67 if getMPISizeWorld() > 10:
68 ne_z = 2*ne_z-1
69 ne_x=int(ceil(ne_z*width_x/depth))
70 #
71 # create receiver array
72 #
73 receiver_line=[2*absorption_zone + i * (rangeRcv//(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
74 #
75 # set source location with tag "source""
76 #
77 src_tags=["source"]
78 src_locations = [ (receiver_line[src_id], depth)]
79 srcloc=(receiver_line[src_id], 0.)
80 #
81 # output
82 #
83 print("%s"%(time.asctime(),))
84 print("ne_x = %s"%(ne_x,))
85 print("ne_z = %s"%(ne_z,))
86 print("width = %s m"%(width_x,))
87 print("depth = %s m"%(depth, ))
88 print("absorption_zone = %s m"%(absorption_zone, ))
89 print("sampling interval = %s ms"%(sampling_interval/U.msec,))
90 print("t_end = %s sec"%(t_end,))
91 print("ricker dominant freqency = %s Hz"%(frq,))
92 print("length of receiver line = %s ms"%(rangeRcv,))
93 print("number of receivers = %s"%(numRcvPerLine,))
94 print("first receiver location = %s m"%(receiver_line[0],))
95 print("last receiver location = %s m"%(receiver_line[-1],))
96 print("source location = %s m"%(src_locations[0][0],))
97 print("source orientation = %s"%(src_dir,))
98 print("matrix lumping = %s"%(lumping,))
99 print("Layer\tV_p\tV_s\teps\tdelta\ttilt\trho")
100 for i in range(len(layers)):
101 print("%s\t%s\t%s\t%s\t%s\t%s\t%s"%( layers[i], v_P[i], v_S[i], eps[i], delta[i], tilt[i], rho[i]))
102 #
103 # create domain:
104 #
105 order = 5
106 domain=Rectangle(order, ne_x,ne_z, l0=width_x, l1=depth,
107 diracPoints=src_locations, diracTags=src_tags, d0=getMPISizeWorld())
108 #
109 # create the wavelet:
110 #
111 wl=Ricker(frq)
112 #
113 #======================================================================
114 #
115 # set
116 #
117 z=ReducedFunction(domain).getX()[1]
118 z_top=0
119 V_P=0
120 V_S=0
121 Delta=0
122 Eps=0
123 Tilt=0
124 Rho=0
125 z_top=depth
126
127 for l in range(len(layers)):
128 m=whereNonPositive(z-z_top)*wherePositive(z-(z_top-layers[l]))
129 V_P = V_P * (1-m) + v_P[l] * m
130 V_S = V_S * (1-m) + v_S[l] * m
131 Delta = Delta * (1-m) + delta[l]* m
132 Eps = Eps * (1-m) + eps[l] * m
133 Tilt = Tilt * (1-m) + tilt[l] * m
134 Rho = Rho * (1-m) + rho[l] * m
135 z_top-=layers[l]
136
137 sw=TTIWave(domain, V_P, V_S, wl, src_tags[0], source_vector = src_dir,
138 eps=Eps, delta=Delta, rho=Rho, theta=Tilt,
139 absorption_zone=absorption_zone, absorption_cut=1e-2, lumping=lumping)
140
141 srclog=Locator(domain, [ (r , depth) for r in receiver_line ] )
142 grploc=[ (x[0], 0.) for x in srclog.getX() ]
143
144 tracer_x=SimpleSEGYWriter(receiver_group=grploc, source=srcloc, sampling_interval=sampling_interval, text='x-displacement')
145 tracer_z=SimpleSEGYWriter(receiver_group=grploc, source=srcloc, sampling_interval=sampling_interval, text='z-displacement')
146
147 if not tracer_x.obspy_available():
148 print("\nWARNING: obspy not available, SEGY files will not be written\n")
149 elif getMPISizeWorld() > 1:
150 print("\nWARNING: SEGY files cannot be written with multiple processes\n")
151
152 t=0.
153 mkDir('output')
154 n=0
155 k_out=0
156 print("calculation starts @ %s"%(time.asctime(),))
157 while t < t_end:
158 t,u = sw.update(t+sampling_interval)
159 tracer_x.addRecord(srclog(u[0]))
160 tracer_z.addRecord(srclog(u[1]))
161 print("t=%s, src=%s: \t %s \t %s \t %s"%(t, wl.getValue(t),srclog(u[1])[0], srclog(u[1])[src_id], srclog(u[1])[-1]))
162 if not n_out is None and n%n_out == 0:
163 print("time step %s written to file %s"%(n_out, "output/u_%d.silo"%(k_out,)))
164 saveSilo("output/u_%d.silo"%(k_out,), u=u)
165 k_out+=1
166 n+=1
167 if tracer_x.obspy_available() and getMPISizeWorld() == 1:
168 tracer_x.write('output/lineX.sgy')
169 tracer_z.write('output/lineZ.sgy')
170 print("calculation completed @ %s"%(time.asctime(),))
171
172 else: # no speckley
173 print("The Speckley module is not available")
174

  ViewVC Help
Powered by ViewVC 1.1.26