/[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 4622 - (show annotations)
Fri Jan 17 04:55:41 2014 UTC (5 years, 3 months ago) by sshaw
File MIME type: text/x-python
File size: 5457 byte(s)
Added dirac support to ripley, added interface for custom assemblers for ripleydomains (also added the custom assembler for 2D VTI waves), changed synthetic_VTI example to use the new, faster custom assembler

1 ##############################################################################
2 #
3 # Copyright (c) 2003-2013 by University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development since 2012 by School of Earth Sciences
11 #
12 ##############################################################################
13 from __future__ import print_function
14
15 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
16 http://www.uq.edu.au
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 from esys.escript import *
23 from esys.escript import unitsSI as U
24 from esys.escript.pdetools import Locator
25 from esys.ripley import Brick, Rectangle
26 from esys.weipa import saveSilo
27 from esys.downunder import Ricker, VTIWave, SimpleSEGYWriter
28 from math import ceil
29
30
31 DIM=2 # spatial dimension
32
33 depth=1*U.km # depth
34 v_p_top=1.5*U.km/U.sec
35 v_p_bottom=3*U.km/U.sec
36 absorption_zone=300*U.m
37 ne_z=500.
38
39 reflector_at=0.5*depth
40
41
42 t_end=0.4*U.sec
43 frq=8.*U.Hz
44 sampling_interval=4*U.msec
45 numRcvPerLine=101
46 rangeRcv=800*U.m
47
48 # location of source in crossing array lines with in 0..numRcvInLine one needs to be None
49 srcEW=numRcvPerLine/2
50 srcNS=None
51
52 # dommain dimension
53 width_x=rangeRcv + 4*absorption_zone
54 width_y=width_x
55 #
56 # create array
57 #
58 receiver_line=[2*absorption_zone + i * (rangeRcv/(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
59 #
60 # set source location with tag "source""
61 #
62 src_tags=["source"]
63 if DIM == 2:
64 src_locations = [ (receiver_line[srcEW], depth)]
65 src_loc_2D=(receiver_line[srcEW], 0.)
66 else:
67 if srcEW:
68 srcNS=numRcvPerLine/2
69 elif srcNS:
70 srcEW=numRcvPerLine/2
71 else:
72 raise ValueError("on of the variables srcEW or srcNS must be None!")
73 src_locations = [ (receiver_line[srcEW], receiver_line[srcNS], depth)]
74 src_loc_2D=(receiver_line[srcEW], receiver_line[srcNS])
75 #
76 # create sensor arrays:
77 #
78 # East-west line of receiver
79 rcvEW_locations=[]
80 rgEW=[]
81 mid_point=receiver_line[len(receiver_line)/2]
82
83 for ix in range(len(receiver_line)):
84 if DIM == 2:
85 rcvEW_locations.append((receiver_line[ix], depth))
86 rgEW.append( ( receiver_line[ix], 0.) )
87 else:
88 rcvEW_locations.append((receiver_line[ix], mid_point, depth))
89 rgEW.append( ( receiver_line[ix], mid_point) )
90 # North-south line of receiver
91 if DIM == 3:
92 rcvNS_locations=[]
93 rgNS=[]
94
95 for iy in range(len(receiver_line)):
96 rcvNS_locations.append((mid_point, receiver_line[iy], depth))
97 rgNS.append( ( mid_point, receiver_line[iy]) )
98 #
99 # create domain:
100 #
101 if DIM == 2:
102 domain=Rectangle(ceil(ne_z*width_x/depth),ne_z,l0=width_x,l1=depth,
103 diracPoints=src_locations, diracTags=src_tags)
104 else:
105 domain=Brick(ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),ne_z,l0=width_x,l1=width_y,l2=depth,
106 diracPoints=src_locations, diracTags=src_tags)
107 wl=Ricker(frq)
108
109 #======================================================================
110 # m=whereNegative(Function(domain).getX()[DIM-1]-reflector_at)
111 # v_p=v_p_bottom*m+v_p_top*(1-m)
112 v_p=2*U.km/U.sec
113 v_s=0.9*U.km/U.sec
114 vareps=0.1*0
115 gamma=0.15*0
116 delta=0.05*0
117 rho=2000*U.kg/U.m**3
118 src_dir=[0,0,1]
119
120 sw=VTIWave(domain, v_p, v_s, wl, src_tags[0], source_vector = src_dir
121 eps=vareps, gamma=gamma, delta=delta, rho=rho,
122 absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True)
123
124 locEW=Locator(domain,rcvEW_locations)
125 tracerEW_x=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='x-displacement - east-west line')
126 tracerEW_z=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='z-displacement - east-west line')
127 if DIM==3:
128 locNS=Locator(domain,rcvNS_locations)
129 tracerEW_y=SimpleSEGYWriter(receiver_group=rgEW, source=src_loc_2D, sampling_interval=sampling_interval, text='x-displacement - east-west line')
130 tracerNS_x=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D, sampling_interval=sampling_interval, text='x-displacement - north-south line')
131 tracerNS_y=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D, sampling_interval=sampling_interval, text='y-displacement - north-south line')
132 tracerNS_z=SimpleSEGYWriter(receiver_group=rgNS, source=src_loc_2D, sampling_interval=sampling_interval, text='z-displacement - north-south line')
133
134
135 t=0.
136 mkDir('tmp')
137 n=0
138 while t < t_end:
139 t,u = sw.update(t+sampling_interval)
140 tracerEW_x.addRecord(locEW(u[0]))
141 tracerEW_z.addRecord(locEW(u[DIM-1]))
142 if DIM==3:
143 tracerEW_y.addRecord(locEW(u[1]))
144 tracerNS_x.addRecord(locNS(u[0]))
145 tracerNS_y.addRecord(locNS(u[1]))
146 tracerNS_z.addRecord(locNS(u[2]))
147 print(t, locEW(u[DIM-1])[len(rgEW)/2-4:len(rgEW)/2+1], wl.getValue(t))
148 #if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n/5,), u=u)
149 saveSilo("tmp/u_%d.silo"%(n,), u=u)
150 n+=1
151 tracerEW_x.write('lineEW_x.sgy')
152 tracerEW_z.write('lineEW_z.sgy')
153 if DIM == 3:
154 tracerEW_y.write('lineEW_y.sgy')
155 tracerNS_x.write('lineNS_x.sgy')
156 tracerNS_y.write('lineNS_y.sgy')
157 tracerNS_z.write('lineNS_z.sgy')

  ViewVC Help
Powered by ViewVC 1.1.26