/[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

trunk/esys2/doc/user/examples/wave.py revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC temp_trunk_copy/doc/examples/wave.py revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC
# Line 1  Line 1 
1  # $Id$  
2    # You can shorten the execution time by reducing variable tend from 60 to 0.5
3    
4  from esys.escript import *  from esys.escript import *
5  from esys.linearPDEs import LinearPDE  from esys.escript.pdetools import Locator
6    from esys.escript.linearPDEs import LinearPDE
7  from esys.finley import Brick  from esys.finley import Brick
8  from numarray import identity  from numarray import identity,zeros,ones
9  ne=10           # number of cells in x_0-direction  
10  depth=10000.   # length in x_0-direction  if not os.path.isdir("data"):
11  width=100000.  # length in x_1 and x_2 direction     print "\nCreating subdirectory 'data'\n"
12       os.mkdir("data")
13    
14    ne=32          # number of cells in x_0 and x_1 directions
15    width=10000.  # length in x_0 and x_1 directions
16  lam=3.462e9  lam=3.462e9
17  mu=3.462e9  mu=3.462e9
18  rho=1154.  rho=1154.
 tau=10.  
 umax=2.  
19  tend=60  tend=60
20  h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
21  print "time step size = ",h  print "time step size = ",h
22    
23  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
24    
25  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):  def wavePropagation(domain,h,tend,lam,mu,rho,U0):
26     x=domain.getX()     x=domain.getX()
27     # ... open new PDE ...     # ... open new PDE ...
28     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
29     mypde.setLumpingOn()     mypde.setSolverMethod(LinearPDE.LUMPING)
30     kronecker=identity(mypde.getDim())     kronecker=identity(mypde.getDim())
31     mypde.setValue(D=kronecker*rho, \  
32                    q=x[0].whereZero()*kronecker[1,:])     #  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     # ... set initial values ....     # ... set initial values ....
44     n=0     n=0
45     u=Vector(0,ContinuousFunction(domain))     # initial value of displacement at point source is constant (U0=0.01)
46     u_last=Vector(0,ContinuousFunction(domain))     # 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     t=0     t=0
50    
51       # define the location of the point source
52       L=Locator(domain,numarray.array(xc))
53       # 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     while t<tend:     while t<tend:
66       # ... get current stress ....       # ... get current stress ....
67       g=grad(u)       g=grad(u)
68       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
69       # ... get new acceleration ....       # ... get new acceleration ....
70       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])       mypde.setValue(X=-stress)          
71       a=mypde.getSolution()       a=mypde.getSolution()
72       # ... get new displacement ...       # ... get new displacement ...
73       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
# Line 45  def wavePropagation(domain,h,tend,lam,mu Line 77  def wavePropagation(domain,h,tend,lam,mu
77       t+=h       t+=h
78       n+=1       n+=1
79       print n,"-th time step t ",t       print n,"-th time step t ",t
80       print "a=",inf(a),sup(a)       u_pc=L.getValue(u)
81       print "u=",inf(u),sup(u)       print "u at point charge=",u_pc
82       # ... save current acceleration in units of gravity      
83       if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10)       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         displacement = length(u), tensor = stress, Ux = u[0] )
93    
94  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)     u_pc_data.close()
95  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)    
96    mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
97    wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
98    

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

  ViewVC Help
Powered by ViewVC 1.1.26