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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26