/[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 113 by jgs, Mon Feb 28 07:06:33 2005 UTC trunk/doc/examples/wave.py revision 2420 by gross, Thu May 14 02:28:58 2009 UTC
# Line 1  Line 1 
1  # $Id$  
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    # You can shorten the execution time by reducing variable tend from 60 to 0.5
23    
24  from esys.escript import *  from esys.escript import *
25  from esys.linearPDEs import LinearPDE  from esys.escript.pdetools import Locator
26    from esys.escript.linearPDEs import LinearPDE
27  from esys.finley import Brick  from esys.finley import Brick
28  from numarray import identity  from numarray import identity,zeros,ones
29  ne=10           # number of cells in x_0-direction  
30  depth=10000.   # length in x_0-direction  if not os.path.isdir("data"):
31  width=100000.  # length in x_1 and x_2 direction     print "\nCreating subdirectory 'data'\n"
32       os.mkdir("data")
33    
34    ne=32          # number of cells in x_0 and x_1 directions
35    width=10000.  # length in x_0 and x_1 directions
36  lam=3.462e9  lam=3.462e9
37  mu=3.462e9  mu=3.462e9
38  rho=1154.  rho=1154.
 tau=10.  
 umax=2.  
39  tend=60  tend=60
40  h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
41  print "time step size = ",h  print "time step size = ",h
42    
43  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
44    
45  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):  def wavePropagation(domain,h,tend,lam,mu,rho,U0):
46     x=domain.getX()     x=domain.getX()
47     # ... open new PDE ...     # ... open new PDE ...
48     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
49     mypde.setLumpingOn()     mypde.setSolverMethod(LinearPDE.LUMPING)
50     kronecker=identity(mypde.getDim())     kronecker=identity(mypde.getDim())
51     mypde.setValue(D=kronecker*rho, \  
52                    q=x[0].whereZero()*kronecker[1,:])     #  spherical source at middle of bottom face
53    
54       xc=[width/2.,width/2.,0.]
55       # define small radius around point xc
56       # Lsup(x) returns the maximum value of the argument x
57       src_radius = 0.1*Lsup(domain.getSize())
58       print "src_radius = ",src_radius
59    
60       dunit=numarray.array([1.,0.,0.]) # defines direction of point source
61    
62       mypde.setValue(D=kronecker*rho)
63     # ... set initial values ....     # ... set initial values ....
64     n=0     n=0
65     u=Vector(0,ContinuousFunction(domain))     # initial value of displacement at point source is constant (U0=0.01)
66     u_last=Vector(0,ContinuousFunction(domain))     # for first two time steps
67       u=U0*whereNegative(length(x-xc)-src_radius)*dunit
68       u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
69     t=0     t=0
70    
71       # define the location of the point source
72       L=Locator(domain,numarray.array(xc))
73       # find potential at point source
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       # open file to save displacement at point source
82       u_pc_data=FileWriter('./data/U_pc.out')
83       u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
84    
85     while t<tend:     while t<tend:
86       # ... get current stress ....       # ... get current stress ....
87       g=grad(u)       g=grad(u)
88       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
89       # ... get new acceleration ....       # ... get new acceleration ....
90       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])       mypde.setValue(X=-stress)          
91       a=mypde.getSolution()       a=mypde.getSolution()
92       # ... get new displacement ...       # ... get new displacement ...
93       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 97  def wavePropagation(domain,h,tend,lam,mu
97       t+=h       t+=h
98       n+=1       n+=1
99       print n,"-th time step t ",t       print n,"-th time step t ",t
100       print "a=",inf(a),sup(a)       u_pc=L.getValue(u)
101       print "u=",inf(u),sup(u)       print "u at point charge=",u_pc
102       # ... save current acceleration in units of gravity      
103       if n%10==0: (length(a)/9.81).saveDX("u.%i.dx"%(n/10))       u_pc_x=u_pc[0]
104         u_pc_y=u_pc[1]
105  print int(width/depth)       u_pc_z=u_pc[2]
106  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)        
107  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)       # save displacements at point source to file for t > 0
108         u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
109    
110         # ... save current acceleration in units of gravity and displacements
111         if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
112         displacement = length(u), tensor = stress, Ux = u[0] )
113    
114       u_pc_data.close()
115      
116    mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
117    wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
118    

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

  ViewVC Help
Powered by ViewVC 1.1.26