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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3148 - (show annotations)
Fri Sep 3 02:09:47 2010 UTC (11 years ago) by jfenwick
File MIME type: text/x-python
File size: 7654 byte(s)
Another attempt to patch the X issue

1
2 ########################################################
3 #
4 # Copyright (c) 2009-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2009-2010 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
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 # You can shorten the execution time by reducing variable tend from 60 to 0.5
23 # Antony Hallam
24 # Acoustic Wave Equation Simulation
25
26 # Importing all the necessary modules required.
27 from esys.escript import *
28 from esys.finley import Rectangle
29 import sys
30 import os
31 from cblib1 import wavesolver2d
32 # smoothing operator
33 from esys.escript.pdetools import Projector
34 import numpy as np
35 import matplotlib
36 matplotlib.use('agg') #It's just here for automated testing
37
38 import pylab as pl
39 import matplotlib.cm as cm
40
41 # Establish a save path.
42 savepath = "data/wavesolver2d009mpltestABCnolump0_0006"
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.0006
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 # ... open new PDE ...
169 mypde=LinearPDE(domain)
170 print mypde.isUsingLumping()
171 print mypde.getSolverOptions()
172 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
173 mypde.setSymmetryOn()
174 kmat = kronecker(domain)
175 mypde.setValue(D=kmat*rho)
176
177 # define small radius around point xc
178 # Lsup(x) returns the maximum value of the argument x
179 src_radius = 50#2*Lsup(domain.getSize())
180 print "src_radius = ",src_radius
181
182 dunit=numpy.array([0.,1.]) # defines direction of point source
183 #~ dunit=(x-xc)
184 #~ absrc=length(dunit)
185 #~ dunit=dunit/maximum(absrc,1e-10)
186
187 # ... set initial values ....
188 n=0
189 # initial value of displacement at point source is constant (U0=0.01)
190 # for first two time steps
191 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
192 #u=whereNegative(length(x-xc)-src_radius)*dunit
193
194 maxi=0.02
195
196 print u
197 u_m1=u
198 t=0
199
200 #~ u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
201 #~ u_pc_x1 = u_pot[0,0]
202 #~ u_pc_y1 = u_pot[0,1]
203 #~ u_pc_x2 = u_pot[1,0]
204 #~ u_pc_y2 = u_pot[1,1]
205 #~ u_pc_x3 = u_pot[2,0]
206 #~ u_pc_y3 = u_pot[2,1]
207 #~
208 #~ # open file to save displacement at point source
209 #~ u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
210 #~ 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))
211
212 while t<tend:
213 # ... get current stress ....
214 # t=1.
215 ##OLD WAY
216 g=grad(u)
217 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
218 ### ... get new acceleration ....
219 #mypde.setValue(X=-stress)
220 #a=mypde.getSolution()
221 ### ... get new displacement ...
222 #u_p1=2*u-u_m1+h*h*a
223 ###NEW WAY
224 mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
225 u_p1 = mypde.getSolution()
226 # ... shift displacements ....
227 u_m1=u
228 u=u_p1*abc
229 #stress =
230 t+=h
231 n+=1
232 print n,"-th time step t ",t
233 #~ u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
234 #~
235 #~ # print "u at point charge=",u_pc
236 #~ u_pc_x1 = u_pot[0,0]
237 #~ u_pc_y1 = u_pot[0,1]
238 #~ u_pc_x2 = u_pot[1,0]
239 #~ u_pc_y2 = u_pot[1,1]
240 #~ u_pc_x3 = u_pot[2,0]
241 #~ u_pc_y3 = u_pot[2,1]
242
243 # save displacements at point source to file for t > 0
244 #~ 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))
245
246 # ... save current acceleration in units of gravity and displacements
247 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
248 #displacement = length(u), tensor = stress, Ux = u[0] )
249 if output == "vtk":
250 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
251 if output == "mpl":
252 uT=np.array(u.toListOfTuples())
253 uT=np.reshape(uT,(ndx+1,ndy+1,2))
254 uTz=uT[:,:,1]+uT[:,:,0]
255 uTz=np.transpose(uTz)
256 pl.clf()
257 # plot wave
258 uTz[0,0]=maxi
259 uTz[0,1]=-maxi
260 CS = pl.imshow(uTz,cmap=cm.spectral)
261 pl.colorbar()
262 # labels and formatting
263 pl.title("Wave Equation Cookbook Example ABC.")
264 pl.xlabel("Horizontal Displacement (m)")
265 pl.ylabel("Depth (m)")
266 if getMPIRankWorld() == 0: #check for MPI processing
267 pl.savefig(os.path.join(savepath,"ws04mpl%05d.png"%n))
268
269 #~ u_pc_data.close()
270 #~ os.system("mencoder mf://"+savepath+"/*.png -mf type=png:\
271 #~ w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
272 #~ wsmpl.avi")
273
274 #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