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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (hide annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 7 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 jfenwick 3981 ##############################################################################
2 ahallam 2829 #
3 jfenwick 5593 # Copyright (c) 2009-2015 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 ahallam 2829 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 sshaw 5707 from __future__ import division, print_function
16 ahallam 2829
17 jfenwick 5593 __copyright__="""Copyright (c) 2009-2015 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ahallam 2829 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 jfenwick 3148 import matplotlib
28     matplotlib.use('agg') #It's just here for automated testing
29    
30 caltinay 4087 import os
31     import sys
32     import numpy as np
33 ahallam 2829 import pylab as pl
34     import matplotlib.cm as cm
35 caltinay 4087 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 ahallam 2829
42     # Establish a save path.
43     savepath = "data/wavesolver2d008mpltestABC"
44 gross 2907 mkDir(savepath)
45 ahallam 2829
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 jfenwick 3892 print("Time step size= ",h, "Expected number of outputs= ",tend/h)
65 ahallam 2829
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 jfenwick 3892 print("src_radius = ",src_radius)
187 ahallam 2829
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 jfenwick 3892 print(u)
211 ahallam 2829 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 jfenwick 3892 print(n,"-th time step t ",t)
252 ahallam 2829 #~ 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