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

revision 327 by gross, Wed Dec 7 04:32:28 2005 UTC revision 582 by lkettle, Wed Mar 8 05:54:06 2006 UTC
# Line 1  Line 1
# \$Id\$
1  from esys.escript import *  from esys.escript import *
2    from esys.escript.pdetools import Locator
3  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
4  from esys.finley import Brick  from esys.finley import Brick
5  ne=5           # number of cells in x_0-direction  from numarray import identity,zeros,ones
6  depth=10000.   # length in x_0-direction
7  width=100000.  # length in x_1 and x_2 direction  ne=32          # number of cells in x_0 and x_1 directions
8    width=10000.  # length in x_0 and x_1 directions
9  lam=3.462e9  lam=3.462e9
10  mu=3.462e9  mu=3.462e9
11  rho=1154.  rho=1154.
tau=10.
umax=2.
12  tend=60  tend=60
13  h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
14  print "time step size = ",h  print "time step size = ",h
15
16  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
17
18  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):  def wavePropagation(domain,h,tend,lam,mu,rho,U0):
19     x=domain.getX()     x=domain.getX()
20     # ... open new PDE ...     # ... open new PDE ...
21     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
22     mypde.setSolverMethod(mypde.LUMPING)     mypde.setSolverMethod(LinearPDE.LUMPING)
23     mypde.setValue(D=kronecker(mypde.getDim())*rho, \     kronecker=identity(mypde.getDim())
24                    q=whereZero(x[0])*kronecker(mypde.getDim())[1,:])
25       #  spherical source at middle of bottom face
26
27       xc=[width/2.,width/2.,0.]
28       # define small radius around point xc
29       # Lsup(x) returns the maximum value of the argument x
32
33       dunit=numarray.array([1.,0.,0.]) # defines direction of point source
34
35       mypde.setValue(D=kronecker*rho)
36     # ... set initial values ....     # ... set initial values ....
37     n=0     n=0
38     u=Vector(0,ContinuousFunction(domain))     # initial value of displacement at point source is constant (U0=0.01)
39     u_last=Vector(0,ContinuousFunction(domain))     # for first two time steps
42     t=0     t=0
43
44       # define the location of the point source
45       L=Locator(domain,xc)
46       # find potential at point source
47       u_pc=L.getValue(u)
48       print "u at point charge=",u_pc
49
50       u_pc_x = u_pc[0]
51       u_pc_y = u_pc[1]
52       u_pc_z = u_pc[2]
53
54       # open file to save displacement at point source
55       u_pc_data=open('./data/U_pc.out','w')
56       u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
57
58     while t<tend:     while t<tend:
59       # ... get current stress ....       # ... get current stress ....
61       stress=lam*trace(g)*kronecker(mypde.getDim())+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
62       # ... get new acceleration ....       # ... get new acceleration ....
63       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker(mypde.getDim())[1,:])       mypde.setValue(X=-stress)
64       a=mypde.getSolution()       a=mypde.getSolution()
65       # ... get new displacement ...       # ... get new displacement ...
66       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 70  def wavePropagation(domain,h,tend,lam,mu
70       t+=h       t+=h
71       n+=1       n+=1
72       print n,"-th time step t ",t       print n,"-th time step t ",t
73       print "a=",inf(a),sup(a)       L=Locator(domain,xc)
74       print "u=",inf(u),sup(u)       u_pc=L.getValue(u)
75       # ... save current acceleration in units of gravity       print "u at point charge=",u_pc
76       if n%10==0: saveVTK("u.%i.xml"%(n/10),a=length(a)/9.81)
77         u_pc_x=u_pc[0]
78         u_pc_y=u_pc[1]
79         u_pc_z=u_pc[2]
80
81         # save displacements at point source to file for t > 0
82         u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
83
84         # ... save current acceleration in units of gravity and displacements
85         if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
86         displacement = length(u), Ux = u[0] )
87
88  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)     u_pc_data.close()
89  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
90    mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
91    wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
92

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