/[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 327 by gross, Wed Dec 7 04:32:28 2005 UTC revision 582 by lkettle, Wed Mar 8 05:54:06 2006 UTC
# Line 1  Line 1 
 # $Id$  
1  from esys.escript import *  from esys.escript import *
2    from esys.escript.pdetools import Locator
3  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
4  from esys.finley import Brick  from esys.finley import Brick
5  ne=5           # number of cells in x_0-direction  from numarray import identity,zeros,ones
6  depth=10000.   # length in x_0-direction  
7  width=100000.  # length in x_1 and x_2 direction  ne=32          # number of cells in x_0 and x_1 directions
8    width=10000.  # length in x_0 and x_1 directions
9  lam=3.462e9  lam=3.462e9
10  mu=3.462e9  mu=3.462e9
11  rho=1154.  rho=1154.
 tau=10.  
 umax=2.  
12  tend=60  tend=60
13  h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
14  print "time step size = ",h  print "time step size = ",h
15    
16  def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)  U0=0.01 # amplitude of point source
17    
18  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):  def wavePropagation(domain,h,tend,lam,mu,rho,U0):
19     x=domain.getX()     x=domain.getX()
20     # ... open new PDE ...     # ... open new PDE ...
21     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
22     mypde.setSolverMethod(mypde.LUMPING)     mypde.setSolverMethod(LinearPDE.LUMPING)
23     mypde.setValue(D=kronecker(mypde.getDim())*rho, \     kronecker=identity(mypde.getDim())
24                    q=whereZero(x[0])*kronecker(mypde.getDim())[1,:])  
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     # ... set initial values ....     # ... set initial values ....
37     n=0     n=0
38     u=Vector(0,ContinuousFunction(domain))     # initial value of displacement at point source is constant (U0=0.01)
39     u_last=Vector(0,ContinuousFunction(domain))     # 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     t=0     t=0
43    
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     while t<tend:     while t<tend:
59       # ... get current stress ....       # ... get current stress ....
60       g=grad(u)       g=grad(u)
61       stress=lam*trace(g)*kronecker(mypde.getDim())+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
62       # ... get new acceleration ....       # ... get new acceleration ....
63       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker(mypde.getDim())[1,:])       mypde.setValue(X=-stress)          
64       a=mypde.getSolution()       a=mypde.getSolution()
65       # ... get new displacement ...       # ... get new displacement ...
66       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
# Line 43  def wavePropagation(domain,h,tend,lam,mu Line 70  def wavePropagation(domain,h,tend,lam,mu
70       t+=h       t+=h
71       n+=1       n+=1
72       print n,"-th time step t ",t       print n,"-th time step t ",t
73       print "a=",inf(a),sup(a)       L=Locator(domain,xc)
74       print "u=",inf(u),sup(u)       u_pc=L.getValue(u)
75       # ... save current acceleration in units of gravity       print "u at point charge=",u_pc
76       if n%10==0: saveVTK("u.%i.xml"%(n/10),a=length(a)/9.81)      
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    
88  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)     u_pc_data.close()
89  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)    
90    mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
91    wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
92    

Legend:
Removed from v.327  
changed lines
  Added in v.582

  ViewVC Help
Powered by ViewVC 1.1.26