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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 113 - (hide annotations)
Mon Feb 28 07:06:33 2005 UTC (14 years, 8 months ago) by jgs
Original Path: trunk/esys2/doc/user/examples/wave.py
File MIME type: text/x-python
File size: 1646 byte(s)
*** empty log message ***

1 jgs 108 # $Id$
2     from esys.escript import *
3 jgs 110 from esys.linearPDEs import LinearPDE
4     from esys.finley import Brick
5     from numarray import identity
6     ne=10 # number of cells in x_0-direction
7     depth=10000. # length in x_0-direction
8     width=100000. # length in x_1 and x_2 direction
9     lam=3.462e9
10     mu=3.462e9
11     rho=1154.
12     tau=10.
13     umax=2.
14     tend=60
15     h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)
16     print "time step size = ",h
17 jgs 108
18 jgs 110 def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)
19    
20     def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):
21 jgs 108 x=domain.getX()
22     # ... open new PDE ...
23 jgs 110 mypde=LinearPDE(domain)
24     mypde.setLumpingOn()
25     kronecker=identity(mypde.getDim())
26     mypde.setValue(D=kronecker*rho, \
27     q=x[0].whereZero()*kronecker[1,:])
28     # ... set initial values ....
29 jgs 108 n=0
30 jgs 110 u=Vector(0,ContinuousFunction(domain))
31     u_last=Vector(0,ContinuousFunction(domain))
32     t=0
33 jgs 108 while t<tend:
34     # ... get current stress ....
35     g=grad(u)
36 jgs 110 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
37 jgs 108 # ... get new acceleration ....
38 jgs 110 mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])
39     a=mypde.getSolution()
40     # ... get new displacement ...
41     u_new=2*u-u_last+h**2*a
42     # ... shift displacements ....
43     u_last=u
44     u=u_new
45     t+=h
46 jgs 108 n+=1
47 jgs 110 print n,"-th time step t ",t
48     print "a=",inf(a),sup(a)
49     print "u=",inf(u),sup(u)
50     # ... save current acceleration in units of gravity
51 jgs 113 if n%10==0: (length(a)/9.81).saveDX("u.%i.dx"%(n/10))
52 jgs 108
53 jgs 113 print int(width/depth)
54 jgs 110 mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)
55     wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
56 jgs 108

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26