/[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 327 - (hide annotations)
Wed Dec 7 04:32:28 2005 UTC (13 years, 2 months ago) by gross
Original Path: trunk/doc/user/examples/wave.py
File MIME type: text/x-python
File size: 1643 byte(s)
make them run again
1 jgs 108 # $Id$
2     from esys.escript import *
3 gross 327 from esys.escript.linearPDEs import LinearPDE
4 jgs 110 from esys.finley import Brick
5 gross 327 ne=5 # number of cells in x_0-direction
6 jgs 110 depth=10000. # length in x_0-direction
7     width=100000. # length in x_1 and x_2 direction
8     lam=3.462e9
9     mu=3.462e9
10     rho=1154.
11     tau=10.
12     umax=2.
13     tend=60
14     h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)
15     print "time step size = ",h
16 jgs 108
17 jgs 110 def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)
18    
19     def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):
20 jgs 108 x=domain.getX()
21     # ... open new PDE ...
22 jgs 110 mypde=LinearPDE(domain)
23 gross 327 mypde.setSolverMethod(mypde.LUMPING)
24     mypde.setValue(D=kronecker(mypde.getDim())*rho, \
25     q=whereZero(x[0])*kronecker(mypde.getDim())[1,:])
26 jgs 110 # ... set initial values ....
27 jgs 108 n=0
28 jgs 110 u=Vector(0,ContinuousFunction(domain))
29     u_last=Vector(0,ContinuousFunction(domain))
30     t=0
31 jgs 108 while t<tend:
32     # ... get current stress ....
33     g=grad(u)
34 gross 327 stress=lam*trace(g)*kronecker(mypde.getDim())+mu*(g+transpose(g))
35 jgs 108 # ... get new acceleration ....
36 gross 327 mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker(mypde.getDim())[1,:])
37 jgs 110 a=mypde.getSolution()
38     # ... get new displacement ...
39     u_new=2*u-u_last+h**2*a
40     # ... shift displacements ....
41     u_last=u
42     u=u_new
43     t+=h
44 jgs 108 n+=1
45 jgs 110 print n,"-th time step t ",t
46     print "a=",inf(a),sup(a)
47     print "u=",inf(u),sup(u)
48     # ... save current acceleration in units of gravity
49 gross 327 if n%10==0: saveVTK("u.%i.xml"%(n/10),a=length(a)/9.81)
50 jgs 108
51 jgs 110 mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)
52     wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
53 jgs 108

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26