/[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 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 6 months 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 ksteube 1215
2     # You can shorten the execution time by reducing variable tend from 60 to 0.5
3    
4 jgs 108 from esys.escript import *
5 lkettle 582 from esys.escript.pdetools import Locator
6 gross 327 from esys.escript.linearPDEs import LinearPDE
7 jgs 110 from esys.finley import Brick
8 lkettle 582 from numarray import identity,zeros,ones
9    
10 ksteube 1215 if not os.path.isdir("data"):
11     print "\nCreating subdirectory 'data'\n"
12     os.mkdir("data")
13    
14 lkettle 582 ne=32 # number of cells in x_0 and x_1 directions
15     width=10000. # length in x_0 and x_1 directions
16 jgs 110 lam=3.462e9
17     mu=3.462e9
18     rho=1154.
19     tend=60
20 lkettle 582 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
21 jgs 110 print "time step size = ",h
22 jgs 108
23 lkettle 582 U0=0.01 # amplitude of point source
24 jgs 110
25 lkettle 582 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
26 jgs 108 x=domain.getX()
27     # ... open new PDE ...
28 jgs 110 mypde=LinearPDE(domain)
29 lkettle 582 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 jgs 110 # ... set initial values ....
44 jgs 108 n=0
45 lkettle 582 # 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 jgs 110 t=0
50 lkettle 582
51     # define the location of the point source
52 jongui 1058 L=Locator(domain,numarray.array(xc))
53 lkettle 582 # 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 jgs 108 while t<tend:
66     # ... get current stress ....
67     g=grad(u)
68 lkettle 582 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
69 jgs 108 # ... get new acceleration ....
70 lkettle 582 mypde.setValue(X=-stress)
71 jgs 110 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 jgs 108 n+=1
79 jgs 110 print n,"-th time step t ",t
80 lkettle 582 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 jongui 1161 displacement = length(u), tensor = stress, Ux = u[0] )
93 jgs 108
94 lkettle 582 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 jgs 108

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26