/[escript]/trunk/doc/examples/usersguide/wave.py
ViewVC logotype

Contents of /trunk/doc/examples/usersguide/wave.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1384 - (show annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 1 month ago) by phornby
Original Path: temp_trunk_copy/doc/examples/wave.py
File MIME type: text/x-python
File size: 2877 byte(s)
Make a temp copy of the trunk before checking in the windows changes


1
2 # You can shorten the execution time by reducing variable tend from 60 to 0.5
3
4 from esys.escript import *
5 from esys.escript.pdetools import Locator
6 from esys.escript.linearPDEs import LinearPDE
7 from esys.finley import Brick
8 from numarray import identity,zeros,ones
9
10 if not os.path.isdir("data"):
11 print "\nCreating subdirectory 'data'\n"
12 os.mkdir("data")
13
14 ne=32 # number of cells in x_0 and x_1 directions
15 width=10000. # length in x_0 and x_1 directions
16 lam=3.462e9
17 mu=3.462e9
18 rho=1154.
19 tend=60
20 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
21 print "time step size = ",h
22
23 U0=0.01 # amplitude of point source
24
25 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
26 x=domain.getX()
27 # ... open new PDE ...
28 mypde=LinearPDE(domain)
29 mypde.setSolverMethod(LinearPDE.LUMPING)
30 kronecker=identity(mypde.getDim())
31
32 # spherical source at middle of bottom face
33
34 xc=[width/2.,width/2.,0.]
35 # define small radius around point xc
36 # Lsup(x) returns the maximum value of the argument x
37 src_radius = 0.1*Lsup(domain.getSize())
38 print "src_radius = ",src_radius
39
40 dunit=numarray.array([1.,0.,0.]) # defines direction of point source
41
42 mypde.setValue(D=kronecker*rho)
43 # ... set initial values ....
44 n=0
45 # initial value of displacement at point source is constant (U0=0.01)
46 # for first two time steps
47 u=U0*whereNegative(length(x-xc)-src_radius)*dunit
48 u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
49 t=0
50
51 # define the location of the point source
52 L=Locator(domain,numarray.array(xc))
53 # find potential at point source
54 u_pc=L.getValue(u)
55 print "u at point charge=",u_pc
56
57 u_pc_x = u_pc[0]
58 u_pc_y = u_pc[1]
59 u_pc_z = u_pc[2]
60
61 # open file to save displacement at point source
62 u_pc_data=open('./data/U_pc.out','w')
63 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
64
65 while t<tend:
66 # ... get current stress ....
67 g=grad(u)
68 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
69 # ... get new acceleration ....
70 mypde.setValue(X=-stress)
71 a=mypde.getSolution()
72 # ... get new displacement ...
73 u_new=2*u-u_last+h**2*a
74 # ... shift displacements ....
75 u_last=u
76 u=u_new
77 t+=h
78 n+=1
79 print n,"-th time step t ",t
80 u_pc=L.getValue(u)
81 print "u at point charge=",u_pc
82
83 u_pc_x=u_pc[0]
84 u_pc_y=u_pc[1]
85 u_pc_z=u_pc[2]
86
87 # save displacements at point source to file for t > 0
88 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
89
90 # ... save current acceleration in units of gravity and displacements
91 if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
92 displacement = length(u), tensor = stress, Ux = u[0] )
93
94 u_pc_data.close()
95
96 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
97 wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
98

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26