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