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/usersguide/wave.py revision 2580 by gross, Tue Aug 4 08:24:12 2009 UTC
# Line 1  Line 1
1  # \$Id\$
2    ########################################################
3    #
4    # Copyright (c) 2003-2009 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  from esys.escript import *  from esys.escript import *
23  from esys.linearPDEs import LinearPDE  from esys.escript.pdetools import Locator
24    from esys.escript.linearPDEs import LinearPDE
25  from esys.finley import Brick  from esys.finley import Brick
26  from numarray import identity  from numpy import identity,zeros,ones
27  ne=10           # number of cells in x_0-direction  import matplotlib.pyplot as plt
28  depth=10000.   # length in x_0-direction
29  width=100000.  # length in x_1 and x_2 direction
30    ne=32          # number of cells in x_0 and x_1 directions
31    width=10000.  # length in x_0 and x_1 directions
32  lam=3.462e9  lam=3.462e9
33  mu=3.462e9  mu=3.462e9
34  rho=1154.  rho=1154.
35  tau=10.  tend=10. # to ran a full simulation change tend to 60.
36  umax=2.  alpha=0.7
37  tend=60  t0=3.
38  h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)
print "time step size = ",h
39
40  def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)  U0=1. # maximum displacement
41    mkDir("data") # create directory data if it does not exist already.
42
43  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):  def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
44       # lists to collect displacement at point source
45       ts, u_pc0,u_pc1,u_pc2=[], [], [], []
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.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
50     kronecker=identity(mypde.getDim())     kronecker=identity(mypde.getDim())
51     mypde.setValue(D=kronecker*rho, \
52                    q=x[0].whereZero()*kronecker[1,:])     dunit=numpy.array([1.,0.,0.]) # defines direction of point source
53
55     # ... set initial values ....     # ... set initial values ....
56     n=0     n=0
57     u=Vector(0,ContinuousFunction(domain))     # for first two time steps
58     u_last=Vector(0,ContinuousFunction(domain))     u=Vector(0.,Solution(domain))
59       u_last=Vector(0.,Solution(domain))
60     t=0     t=0
61
62       # define the location of the point source
63       L=Locator(domain,xc)
64       # find potential at point source
65       u_pc=L.getValue(u)
66       print "u at point charge=",u_pc
67       ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
68
69     while t<tend:     while t<tend:
70         t+=h
71       # ... get current stress ....       # ... get current stress ....
73       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
74       # ... get new acceleration ....       # ... get new acceleration ....
75       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])       amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2*exp(1./2.-(t-t0)**2/alpha**2)
76         mypde.setValue(X=-stress, r=dunit*amplitude)
77       a=mypde.getSolution()       a=mypde.getSolution()
78       # ... get new displacement ...       # ... get new displacement ...
79       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
80       # ... shift displacements ....       # ... shift displacements ....
81       u_last=u       u_last=u
82       u=u_new       u=u_new
t+=h
83       n+=1       n+=1
84       print n,"-th time step t ",t       print n,"-th time step t ",t
85       print "a=",inf(a),sup(a)       u_pc=L.getValue(u)
86       print "u=",inf(u),sup(u)       print "u at point charge=",u_pc
87       # ... save current acceleration in units of gravity       ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
88       if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10)
89         # ... save current acceleration in units of gravity and displacements
90         if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
91         displacement = length(u), tensor = stress, Ux = u[0] )
92       return ts, u_pc0,u_pc1,u_pc2
93
94  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)  #
95  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)  # create domain:
96    #
97    mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
98    #
99    #  sety time step size:
100    #
101    h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
102    print "time step size = ",h
103    #
104    #  spherical source at middle of bottom face
105    #
106    xc=[width/2.,width/2.,0.]
107    # define small radius around point xc
110    #
111    # run it
112    #
113    ts, u_pc0,u_pc1,u_pc2 = wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)
114    #
115    # create a plot:
116    #
117    if getMPIRankWorld() == 0:
118        plt.title("Displacement at Point Source")
119        plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
120        plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
121        plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
122        plt.xlabel('time')
123        plt.ylabel('displacement')
124        plt.legend()
125        plt.savefig('u_pc.png', format='png')
126    # or save displacement
127    u_pc_data=FileWriter('./data/U_pc.out')
128    for i in xrange(len(ts)) :
129        u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
130    u_pc_data.close()
131

Legend:
 Removed from v.110 changed lines Added in v.2580