/[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 108 by jgs, Thu Jan 27 06:21:59 2005 UTC trunk/doc/user/examples/wave.py 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  import esys.finley  from esys.escript.pdetools import Locator
3  import numarray  from esys.escript.linearPDEs import LinearPDE
4    from esys.finley import Brick
5    from numarray import identity,zeros,ones
6    
7  def wavePropagation(domain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax):  ne=32          # number of cells in x_0 and x_1 directions
8     # ... get characteristic function of impact region:  width=10000.  # length in x_0 and x_1 directions
9    lam=3.462e9
10    mu=3.462e9
11    rho=1154.
12    tend=60
13    h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
14    print "time step size = ",h
15    
16    U0=0.01 # amplitude of point source
17    
18    def wavePropagation(domain,h,tend,lam,mu,rho,U0):
19     x=domain.getX()     x=domain.getX()
    location=x[0].whereZero()*(length(x-xc)-r).whereNegative()  
20     # ... open new PDE ...     # ... open new PDE ...
21     myPDE=LinearPDE(mydomain)     mypde=LinearPDE(domain)
22     myPDE.setLumpingOn()     mypde.setSolverMethod(LinearPDE.LUMPING)
23     myPDE.setValue(D=numarray.identity(myPDE.getDim())*rho,q=location*numarray.identity(myPDE.getDim())[:,0])     kronecker=identity(mypde.getDim())
24     # ... set the initial values :  
25     t=dt     #  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 ....
37     n=0     n=0
38     u=0     # initial value of displacement at point source is constant (U0=0.01)
39     v=0     # for first two time steps
40     a=0     u=U0*whereNegative(length(x-xc)-src_radius)*dunit
41       u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
42       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:
      # ... up-date displacement ....    
      u=u+dt*v+dt**2/2*a  
59       # ... get current stress ....       # ... get current stress ....
60       g=grad(u)       g=grad(u)
61       stress=lame_lambda*trace(g)+lame_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,q=impact_location, \       mypde.setValue(X=-stress)          
64                  r=umax/tau**2*(6*t/tau-9*(t/tau)^4)*exp(-(t/tau)^3)*numarray.identity(myPDE.getDim())[:,0])       a=mypde.getSolution()
65       a_new=myPDE.getSolution()       # ... get new displacement ...
66       # ... update velocity ...       u_new=2*u-u_last+h**2*a
67       v=v+h/2*(a+a_new)       # ... shift displacements ....
68       # ... next time step ...       u_last=u
69       a=a_new       u=u_new
70       t+=dt       t+=h
71       n+=1       n+=1
72       # ... save current displacement:       print n,"-th time step t ",t
73       if n%10: u.saveDX("u.%i.dx"%n)       L=Locator(domain,xc)
74         u_pc=L.getValue(u)
75         print "u at point charge=",u_pc
76        
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  ne=6     u_pc_data.close()
89  lame_lambda=3.462e9    
90  lame_mu=3.462e9  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
91  rho=1154.  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
 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)  
92    

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

  ViewVC Help
Powered by ViewVC 1.1.26