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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2  from esys.escript import *  from esys.escript import *
3  import esys.finley  from esys.linearPDEs import LinearPDE
4  import numarray  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    
18    def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)
19    
20  def wavePropagation(domain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax):  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):
    # ... get characteristic function of impact region:  
21     x=domain.getX()     x=domain.getX()
    location=x[0].whereZero()*(length(x-xc)-r).whereNegative()  
22     # ... open new PDE ...     # ... open new PDE ...
23     myPDE=LinearPDE(mydomain)     mypde=LinearPDE(domain)
24     myPDE.setLumpingOn()     mypde.setLumpingOn()
25     myPDE.setValue(D=numarray.identity(myPDE.getDim())*rho,q=location*numarray.identity(myPDE.getDim())[:,0])     kronecker=identity(mypde.getDim())
26     # ... set the initial values :     mypde.setValue(D=kronecker*rho, \
27     t=dt                    q=x[0].whereZero()*kronecker[1,:])
28       # ... set initial values ....
29     n=0     n=0
30     u=0     u=Vector(0,ContinuousFunction(domain))
31     v=0     u_last=Vector(0,ContinuousFunction(domain))
32     a=0     t=0
33     while t<tend:     while t<tend:
      # ... up-date displacement ....    
      u=u+dt*v+dt**2/2*a  
34       # ... get current stress ....       # ... get current stress ....
35       g=grad(u)       g=grad(u)
36       stress=lame_lambda*trace(g)+lame_mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
37       # ... get new acceleration ....       # ... get new acceleration ....
38       myPDE.setValue(X=stress,q=impact_location, \       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])
39                  r=umax/tau**2*(6*t/tau-9*(t/tau)^4)*exp(-(t/tau)^3)*numarray.identity(myPDE.getDim())[:,0])       a=mypde.getSolution()
40       a_new=myPDE.getSolution()       # ... get new displacement ...
41       # ... update velocity ...       u_new=2*u-u_last+h**2*a
42       v=v+h/2*(a+a_new)       # ... shift displacements ....
43       # ... next time step ...       u_last=u
44       a=a_new       u=u_new
45       t+=dt       t+=h
46       n+=1       n+=1
47       # ... save current displacement:       print n,"-th time step t ",t
48       if n%10: u.saveDX("u.%i.dx"%n)       print "a=",inf(a),sup(a)
49         print "u=",inf(u),sup(u)
50         # ... save current acceleration in units of gravity
51         if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10)
52    
53  ne=6  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)
54  lame_lambda=3.462e9  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
 lame_mu=3.462e9  
 rho=1154.  
 tau=2.  
 umax=15.  
 xc=[0,1000,1000]  
 r=1.  
 tend=10.  
 dt=1./5.*sqrt(rho/(lame_lambda+2*lame_mu))(20000./ne)  
 print "step size = ",dt  
 mydomain=finely.Brick(ne,10*ne,10*ne,l0=20000,l1=200000,l2=200000)  
 wavePropagation(mydomain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax)  
55    

Legend:
Removed from v.108  
changed lines
  Added in v.110

  ViewVC Help
Powered by ViewVC 1.1.26