/[escript]/trunk/doc/examples/cookbook/wavesolver2d004.py
ViewVC logotype

Contents of /trunk/doc/examples/cookbook/wavesolver2d004.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 6 months ago) by sshaw
File MIME type: text/x-python
File size: 8068 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2009-2015 by The 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 until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15 from __future__ import division, print_function
16
17 __copyright__="""Copyright (c) 2009-2015 by The University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 # You can shorten the execution time by reducing variable tend from 60 to 0.5
25
26 # Importing all the necessary modules required.
27 import matplotlib
28 matplotlib.use('agg') #It's just here for automated testing
29
30 import os
31 import sys
32 import numpy as np
33 import pylab as pl
34 import matplotlib.cm as cm
35 from esys.escript import *
36 # smoothing operator
37 from esys.escript.pdetools import Projector
38 from esys.finley import Rectangle
39 from esys.weipa import saveVTK
40 from cblib1 import wavesolver2d
41
42 # Establish a save path.
43 savepath = "data/wavesolver2d008mpltestABC"
44 mkDir(savepath)
45
46
47 #Geometric and material property related variables.
48 mx = 1000. # model lenght
49 my = 1000. # model width
50 ndx = 200 # steps in x direction
51 ndy = 200 # steps in y direction
52
53 xstep=mx/ndx
54 ystep=my/ndy
55
56 lam=3.462e9 #lames constant
57 mu=3.462e9 #bulk modulus
58 rho=1154. #density
59 # Time related variables.
60 tend=0.5 #end time
61 #calculating )the timestep
62 h=(1./5.)*sqrt(rho/(lam+2*mu))*(mx/ndx)
63 #Check to make sure number of time steps is not too large.
64 print("Time step size= ",h, "Expected number of outputs= ",tend/h)
65
66 #uncomment the following lines to give the user a chance to stop
67 #proceeder = raw_input("Is this ok?(y/n)")
68 #Exit if user thinks too many outputs.
69 #if proceeder == "n":
70 # sys.exit()
71
72 U0=0.01 # amplitude of point source
73 # spherical source at middle of bottom face
74
75 xc=[500,500]
76
77 mydomain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
78 #wavesolver2d(mydomain,h,tend,lam,mu,rho,U0,xc,savepath,output="mpl")
79
80
81
82
83 domain=mydomain
84 output="mpl"
85
86
87
88
89
90 from esys.escript.linearPDEs import LinearPDE
91 x=domain.getX()
92
93 ## boundary conditions
94
95 bleft=xstep*50.
96 bright=mx-(xstep*50.)
97 bbot=my-(ystep*50.)
98 btop=ystep*50.
99
100 left=x[0]-bleft
101 right=x[0]-bright
102 bottom=x[1]-bbot
103 top=x[1]-btop
104
105 decay=0.0005
106 fleft=exp(-1.0*(decay*(bleft-x[0]))**2)
107 fright=exp(-1.0*(decay*(x[0]-bright))**2)
108 fbottom=exp(-1.0*(decay*(x[1]-bbot))**2)
109 ftop=exp(-1.0*(decay*(btop-x[1]))**2)
110
111 abcleft=fleft*whereNegative(left)
112 abcright=fright*wherePositive(right)
113 abcbottom=fbottom*wherePositive(bottom)
114 abctop=ftop*whereNegative(top)
115
116 abcleft=abcleft+whereZero(abcleft)
117 abcright=abcright+whereZero(abcright)
118 abcbottom=abcbottom+whereZero(abcbottom)
119 abctop=abctop+whereZero(abctop)
120
121 abc=abcleft*abcright*abcbottom*abctop
122
123 #~ fleftT=fleft.toListOfTuples()
124 #~ fleftT=np.reshape(fleftT,(ndx+1,ndy+1))
125 #~ pl.imshow(fleftT)
126 #~ pl.colorbar()
127 #~ pl.savefig("fleftT.png")
128 #~
129 #~ frightT=fright.toListOfTuples()
130 #~ frightT=np.reshape(frightT,(ndx+1,ndy+1))
131 #~ pl.clf()
132 #~ pl.imshow(frightT)
133 #~ pl.colorbar()
134 #~ pl.savefig("frightT.png")
135 #~
136 #~ fbottomT=fbottom.toListOfTuples()
137 #~ fbottomT=np.reshape(fbottomT,(ndx+1,ndy+1))
138 #~ pl.clf()
139 #~ pl.imshow(fbottomT)
140 #~ pl.colorbar()
141 #~ pl.savefig("fbottomT.png")
142 #~
143 #~ #tester=fright*wherePositive(right)
144 #~ tester=fleft*whereNegative(left)
145 #~ tester=tester.toListOfTuples()
146 #~ tester=np.reshape(tester,(ndx+1,ndy+1))
147 #~ pl.clf()
148 #~ pl.imshow(tester)
149 #~ pl.colorbar()
150 #~ pl.savefig("tester1.png")
151 #~
152 #~ tester=fright*wherePositive(right)
153 #~ tester=tester.toListOfTuples()
154 #~ tester=np.reshape(tester,(ndx+1,ndy+1))
155 #~ pl.clf()
156 #~ pl.imshow(tester)
157 #~ pl.colorbar()
158 #~ pl.savefig("tester2.png")
159 #~
160 #~ tester=fbottom*wherePositive(bottom)
161 #~ tester=tester.toListOfTuples()
162 #~ tester=np.reshape(tester,(ndx+1,ndy+1))
163 #~ pl.clf()
164 #~ pl.imshow(tester)
165 #~ pl.colorbar()
166 #~ pl.savefig("tester3.png")
167
168
169 abcT=abc.toListOfTuples()
170 abcT=np.reshape(abcT,(ndx+1,ndy+1))
171 pl.clf()
172 pl.imshow(abcT)
173 pl.colorbar()
174 pl.savefig("abc.png")
175
176 # ... open new PDE ...
177 mypde=LinearPDE(domain)
178 #mypde.setSolverMethod(LinearPDE.LUMPING)
179 mypde.setSymmetryOn()
180 kmat = kronecker(domain)
181 mypde.setValue(D=kmat*rho)
182
183 # define small radius around point xc
184 # Lsup(x) returns the maximum value of the argument x
185 src_radius = 50#2*Lsup(domain.getSize())
186 print("src_radius = ",src_radius)
187
188 #dunit=numpy.array([0.,1.]) # defines direction of point source
189 dunit=(x-xc)
190 absrc=length(dunit)
191 dunit=dunit/maximum(absrc,1e-10)
192
193 # ... set initial values ....
194 n=0
195 # initial value of displacement at point source is constant (U0=0.01)
196 # for first two time steps
197 u=x*0.
198 #u=abs(U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit)
199 #u=whereNegative(length(x-xc)-src_radius)*dunit
200
201 #maxi=max(u.toListOfTuples())
202
203
204 #~ srcf
205 #~
206 #~ x2= np.linspace(0,0.1,333)
207 #~
208 #~ y2=np.exp(-50.*x2)*np.sin(40*3.14157*x2)
209
210 print(u)
211 u_m1=u
212 t=0
213
214 #~ u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
215 #~ u_pc_x1 = u_pot[0,0]
216 #~ u_pc_y1 = u_pot[0,1]
217 #~ u_pc_x2 = u_pot[1,0]
218 #~ u_pc_y2 = u_pot[1,1]
219 #~ u_pc_x3 = u_pot[2,0]
220 #~ u_pc_y3 = u_pot[2,1]
221 #~
222 #~ # open file to save displacement at point source
223 #~ u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
224 #~ u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
225
226 while t<tend:
227 # ... get current stress ....
228 # t=1.
229
230 ## TIMED SOURCE
231 srcf=np.exp(-50.*t)*np.sin(40*3.14157*t)
232 setValueOfDataPoint(srcf,xc,u)
233 ##OLD WAY
234 g=grad(u)
235 break
236 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
237 ### ... get new acceleration ....
238 #mypde.setValue(X=-stress)
239 #a=mypde.getSolution()
240 ### ... get new displacement ...
241 #u_p1=2*u-u_m1+h*h*a
242 ###NEW WAY
243 mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
244 u_p1 = mypde.getSolution()
245 # ... shift displacements ....
246 u_m1=u
247 u=u_p1*abc
248 #stress =
249 t+=h
250 n+=1
251 print(n,"-th time step t ",t)
252 #~ u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
253 #~
254 #~ # print "u at point charge=",u_pc
255 #~ u_pc_x1 = u_pot[0,0]
256 #~ u_pc_y1 = u_pot[0,1]
257 #~ u_pc_x2 = u_pot[1,0]
258 #~ u_pc_y2 = u_pot[1,1]
259 #~ u_pc_x3 = u_pot[2,0]
260 #~ u_pc_y3 = u_pot[2,1]
261
262 # save displacements at point source to file for t > 0
263 #~ u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
264
265 # ... save current acceleration in units of gravity and displacements
266 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
267 #displacement = length(u), tensor = stress, Ux = u[0] )
268 if output == "vtk":
269 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
270 if output == "mpl":
271 uT=np.array(u.toListOfTuples())
272 uT=np.reshape(uT,(ndx+1,ndy+1,2))
273 uTz=uT[:,:,1]+uT[:,:,0]
274 uTz=np.transpose(uTz)
275 pl.clf()
276 # plot wave
277 uTz[0,0]=maxi
278 uTz[0,1]=-maxi
279 CS = pl.imshow(uTz,cmap=cm.spectral)
280 pl.colorbar()
281 # labels and formatting
282 pl.title("Wave Equation Cookbook Example ABC.")
283 pl.xlabel("Horizontal Displacement (m)")
284 pl.ylabel("Depth (m)")
285 if getMPIRankWorld() == 0: #check for MPI processing
286 pl.savefig(os.path.join(savepath,"ws04mpl%05d.png"%n))
287
288 #~ u_pc_data.close()
289 #~ os.system("mencoder mf://"+savepath+"/*.png -mf type=png:\
290 #~ w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
291 #~ wsmpl.avi")
292
293 #mencoder mf://*.png -mf type=png:\w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o wsmpl.avi

  ViewVC Help
Powered by ViewVC 1.1.26