/[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/doc/examples/wave.py revision 2420 by gross, Thu May 14 02:28:58 2009 UTC trunk/doc/examples/usersguide/wave.py revision 2549 by jfenwick, Mon Jul 20 06:43:47 2009 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 19  __license__="""Licensed under the Open S Line 19  __license__="""Licensed under the Open S
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
 # You can shorten the execution time by reducing variable tend from 60 to 0.5  
   
22  from esys.escript import *  from esys.escript import *
23  from esys.escript.pdetools import Locator  from esys.escript.pdetools import Locator
24  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
25  from esys.finley import Brick  from esys.finley import Brick
26  from numarray import identity,zeros,ones  from numpy import identity,zeros,ones
   
 if not os.path.isdir("data"):  
    print "\nCreating subdirectory 'data'\n"  
    os.mkdir("data")  
27    
28  ne=32          # number of cells in x_0 and x_1 directions  ne=32          # number of cells in x_0 and x_1 directions
29  width=10000.  # length in x_0 and x_1 directions  width=10000.  # length in x_0 and x_1 directions
30  lam=3.462e9  lam=3.462e9
31  mu=3.462e9  mu=3.462e9
32  rho=1154.  rho=1154.
33  tend=60  tend=20. # to ran a full simulation change tend to 60.
34  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)  alpha=0.7
35  print "time step size = ",h  t0=3.
36    
37  U0=0.01 # amplitude of point source  U0=1. # maximum displacement
38    mkDir("data") # create directory data if it does not exist already.
39    
40  def wavePropagation(domain,h,tend,lam,mu,rho,U0):  def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
41     x=domain.getX()     x=domain.getX()
42     # ... open new PDE ...     # ... open new PDE ...
43     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
44     mypde.setSolverMethod(LinearPDE.LUMPING)     mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
45     kronecker=identity(mypde.getDim())     kronecker=identity(mypde.getDim())
46    
47     #  spherical source at middle of bottom face     dunit=numpy.array([1.,0.,0.]) # defines direction of point source
   
    xc=[width/2.,width/2.,0.]  
    # define small radius around point xc  
    # Lsup(x) returns the maximum value of the argument x  
    src_radius = 0.1*Lsup(domain.getSize())  
    print "src_radius = ",src_radius  
   
    dunit=numarray.array([1.,0.,0.]) # defines direction of point source  
48    
49     mypde.setValue(D=kronecker*rho)     mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
50     # ... set initial values ....     # ... set initial values ....
51     n=0     n=0
    # initial value of displacement at point source is constant (U0=0.01)  
52     # for first two time steps     # for first two time steps
53     u=U0*whereNegative(length(x-xc)-src_radius)*dunit     u=Vector(0.,Solution(domain))
54     u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit     u_last=Vector(0.,Solution(domain))
55     t=0     t=0
56    
57     # define the location of the point source     # define the location of the point source
58     L=Locator(domain,numarray.array(xc))     L=Locator(domain,xc)
59     # find potential at point source     # find potential at point source
60     u_pc=L.getValue(u)     u_pc=L.getValue(u)
61     print "u at point charge=",u_pc     print "u at point charge=",u_pc
     
    u_pc_x = u_pc[0]  
    u_pc_y = u_pc[1]  
    u_pc_z = u_pc[2]  
   
62     # open file to save displacement at point source     # open file to save displacement at point source
63     u_pc_data=FileWriter('./data/U_pc.out')     u_pc_data=FileWriter('./data/U_pc.out')
64     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))     u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
65    
66     while t<tend:     while t<tend:
67         t+=h
68       # ... get current stress ....       # ... get current stress ....
69       g=grad(u)       g=grad(u)
70       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
71       # ... get new acceleration ....       # ... get new acceleration ....
72       mypde.setValue(X=-stress)                 amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2*exp(1./2.-(t-t0)**2/alpha**2)
73         mypde.setValue(X=-stress, r=dunit*amplitude)
74       a=mypde.getSolution()       a=mypde.getSolution()
75       # ... get new displacement ...       # ... get new displacement ...
76       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
77       # ... shift displacements ....       # ... shift displacements ....
78       u_last=u       u_last=u
79       u=u_new       u=u_new
      t+=h  
80       n+=1       n+=1
81       print n,"-th time step t ",t       print n,"-th time step t ",t
82       u_pc=L.getValue(u)       u_pc=L.getValue(u)
83       print "u at point charge=",u_pc       print "u at point charge=",u_pc
       
      u_pc_x=u_pc[0]  
      u_pc_y=u_pc[1]  
      u_pc_z=u_pc[2]  
         
84       # save displacements at point source to file for t > 0       # save displacements at point source to file for t > 0
85       u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))       u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
86    
87       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
88       if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,       if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
# Line 113  def wavePropagation(domain,h,tend,lam,mu Line 90  def wavePropagation(domain,h,tend,lam,mu
90    
91     u_pc_data.close()     u_pc_data.close()
92        
93  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
94  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)  h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
95    print "time step size = ",h
96    #  spherical source at middle of bottom face
97    xc=[width/2.,width/2.,0.]
98    # define small radius around point xc
99    src_radius = 0.03*width
100    print "src_radius = ",src_radius
101    wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)
102    

Legend:
Removed from v.2420  
changed lines
  Added in v.2549

  ViewVC Help
Powered by ViewVC 1.1.26