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

Annotation of /trunk/doc/examples/inversion/synthetic_VTI.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6812 - (hide annotations)
Fri May 3 00:36:47 2019 UTC (7 months, 1 week ago) by aellery
File MIME type: text/x-python
File size: 6677 byte(s)
Another silo fix


1 gross 4548 ##############################################################################
2     #
3 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
4 gross 4548 # http://www.uq.edu.au
5     #
6     # Primary Business: Queensland, Australia
7 jfenwick 6112 # Licensed under the Apache License, version 2.0
8     # http://www.apache.org/licenses/LICENSE-2.0
9 gross 4548 #
10 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
11     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
12 gross 4548 #
13     ##############################################################################
14 sshaw 5666 from __future__ import print_function, division
15 gross 4548
16 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
17 gross 4548 http://www.uq.edu.au
18     Primary Business: Queensland, Australia"""
19 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
20     http://www.apache.org/licenses/LICENSE-2.0"""
21 gross 4548 __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 gross 4562 from esys.downunder import Ricker, VTIWave, SimpleSEGYWriter
28 gross 4548 from math import ceil
29    
30 caltinay 6119 try:
31     from esys.speckley import Brick, Rectangle
32     HAVE_SPECKLEY=True
33     except ImportError:
34     HAVE_SPECKLEY=False
35 gross 4548
36 caltinay 6119 if HAVE_SPECKLEY:
37     DIM=2 # spatial dimension
38 gross 4548
39 caltinay 6119 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 gross 4553
45 caltinay 6119 reflector_at=0.5*depth
46 gross 4548
47    
48 caltinay 6119 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 gross 4548
54 caltinay 6119 # location of source in crossing array lines with in 0..numRcvInLine one needs to be None
55     srcEW=numRcvPerLine//2
56     srcNS=None
57 gross 4548
58 caltinay 6119 # 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 gross 4553
89 caltinay 6119 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 gross 4553
101 caltinay 6119 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 gross 4548
117 caltinay 6119 #======================================================================
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 gross 4548
128 caltinay 6119 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 gross 4562
132 caltinay 6119 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 gross 4548
154 caltinay 6119 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 aellery 6812 saveSilo("tmp/u_%d.silo"%(n,), u=u, cycle=n, time=t)
169 caltinay 6119 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