# Diff of /trunk/doc/examples/usersguide/wave.py

trunk/doc/user/examples/wave.py revision 327 by gross, Wed Dec 7 04:32:28 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
11    #
12    ########################################################
13
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
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.escript.pdetools import Locator
26  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
27  from esys.finley import Brick  from esys.finley import Brick
28  ne=5           # number of cells in x_0-direction  from numarray import identity,zeros,ones
29  depth=10000.   # length in x_0-direction
30  width=100000.  # length in x_1 and x_2 direction  if not os.path.isdir("data"):
31       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.setSolverMethod(mypde.LUMPING)     mypde.setSolverMethod(LinearPDE.LUMPING)
50     mypde.setValue(D=kronecker(mypde.getDim())*rho, \     kronecker=identity(mypde.getDim())
51                    q=whereZero(x[0])*kronecker(mypde.getDim())[1,:])
52       #  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
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
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 ....
88       stress=lam*trace(g)*kronecker(mypde.getDim())+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(mypde.getDim())[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 43  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: saveVTK("u.%i.xml"%(n/10),a=length(a)/9.81)       u_pc_x=u_pc[0]
104         u_pc_y=u_pc[1]
105         u_pc_z=u_pc[2]
106
107         # 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  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)     u_pc_data.close()
115  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
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.327 changed lines Added in v.2420