/[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/user/examples/wave.py revision 582 by lkettle, Wed Mar 8 05:54:06 2006 UTC trunk/doc/examples/usersguide/wave.py revision 2502 by gross, Tue Jun 30 05:49:22 2009 UTC
# Line 1  Line 1 
1    
2    ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="https://launchpad.net/escript-finley"
21    
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
27    
28    if not os.path.isdir("data"):
29       print "\nCreating subdirectory 'data'\n"
30       if getMPIRankWorld()==0:  os.mkdir("data")
31    
32  ne=32          # number of cells in x_0 and x_1 directions  ne=32          # number of cells in x_0 and x_1 directions
33  width=10000.  # length in x_0 and x_1 directions  width=10000.  # length in x_0 and x_1 directions
34  lam=3.462e9  lam=3.462e9
35  mu=3.462e9  mu=3.462e9
36  rho=1154.  rho=1154.
37  tend=60  tend=0.05     # to ran a full simulation change tend to 60.
38  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)  alpha=0.3
 print "time step size = ",h  
39    
40  U0=0.01 # amplitude of point source  U0=0.01 # amplitude of point source
41    
# Line 19  def wavePropagation(domain,h,tend,lam,mu Line 43  def wavePropagation(domain,h,tend,lam,mu
43     x=domain.getX()     x=domain.getX()
44     # ... open new PDE ...     # ... open new PDE ...
45     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
46     mypde.setSolverMethod(LinearPDE.LUMPING)     mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
47     kronecker=identity(mypde.getDim())     kronecker=identity(mypde.getDim())
48    
49     #  spherical source at middle of bottom face     #  spherical source at middle of bottom face
50    
51     xc=[width/2.,width/2.,0.]     xc=[width/2.,width/2.,0.]
52     # define small radius around point xc     # define small radius around point xc
53     # Lsup(x) returns the maximum value of the argument x     src_radius = 0.03*width
    src_radius = 0.1*Lsup(domain.getSize())  
54     print "src_radius = ",src_radius     print "src_radius = ",src_radius
55    
56     dunit=numarray.array([1.,0.,0.]) # defines direction of point source     dunit=numpy.array([1.,0.,0.]) # defines direction of point source
57    
58     mypde.setValue(D=kronecker*rho)     mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
59     # ... set initial values ....     # ... set initial values ....
60     n=0     n=0
61     # initial value of displacement at point source is constant (U0=0.01)     # initial value of displacement at point source is constant (U0=0.01)
62     # for first two time steps     # for first two time steps
63     u=U0*whereNegative(length(x-xc)-src_radius)*dunit     u=Vector(0.,Solution(domain))
64     u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit     u_last=Vector(0.,Solution(domain))
65     t=0     t=0
66    
67     # define the location of the point source     # define the location of the point source
# Line 46  def wavePropagation(domain,h,tend,lam,mu Line 69  def wavePropagation(domain,h,tend,lam,mu
69     # find potential at point source     # find potential at point source
70     u_pc=L.getValue(u)     u_pc=L.getValue(u)
71     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]  
   
72     # open file to save displacement at point source     # open file to save displacement at point source
73     u_pc_data=open('./data/U_pc.out','w')     u_pc_data=FileWriter('./data/U_pc.out')
74     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]))
75    
76     while t<tend:     while t<tend:
77         t+=h
78       # ... get current stress ....       # ... get current stress ....
79       g=grad(u)       g=grad(u)
80       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
81       # ... get new acceleration ....       # ... get new acceleration ....
82       mypde.setValue(X=-stress)                 amplitude=U0*2*exp(1)/alpha**2*(1-5*(t/alpha)**2+2*(t/alpha)**4)*exp(-(t/alpha)**2)
83         mypde.setValue(X=-stress, r=dunit*amplitude)
84       a=mypde.getSolution()       a=mypde.getSolution()
85       # ... get new displacement ...       # ... get new displacement ...
86       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
87       # ... shift displacements ....       # ... shift displacements ....
88       u_last=u       u_last=u
89       u=u_new       u=u_new
      t+=h  
90       n+=1       n+=1
91       print n,"-th time step t ",t       print n,"-th time step t ",t
      L=Locator(domain,xc)  
92       u_pc=L.getValue(u)       u_pc=L.getValue(u)
93       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]  
         
94       # save displacements at point source to file for t > 0       # save displacements at point source to file for t > 0
95       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]))
96    
97       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
98       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,
99       displacement = length(u), Ux = u[0] )       displacement = length(u), tensor = stress, Ux = u[0] )
100    
101     u_pc_data.close()     u_pc_data.close()
102        
103  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)
104    h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
105    print "time step size = ",h
106  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
107    

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

  ViewVC Help
Powered by ViewVC 1.1.26