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