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

trunk/esys2/doc/user/examples/wave.py revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC trunk/doc/examples/wave.py revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 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
20    __url__="http://www.uq.edu.au/esscc/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
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=open('./data/U_pc.out','w')
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+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("/tmp/res/u.%i.dx"%(n/10)       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.110 changed lines Added in v.1811