/[escript]/trunk/doc/examples/usersguide/wave.py
ViewVC logotype

Annotation of /trunk/doc/examples/usersguide/wave.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years ago) by jfenwick
File MIME type: text/x-python
File size: 4274 byte(s)
Don't panic.
Updating copyright stamps

1 ksteube 1215
2 ksteube 1811 ########################################################
3     #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 # 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 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1811 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1811
22 jgs 108 from esys.escript import *
23 lkettle 582 from esys.escript.pdetools import Locator
24 gross 327 from esys.escript.linearPDEs import LinearPDE
25 jgs 110 from esys.finley import Brick
26 jfenwick 2455 from numpy import identity,zeros,ones
27 jfenwick 2600 import matplotlib
28    
29     matplotlib.use('agg') #For interactive use, you can comment out this line
30     #It's just here to make testing easier
31    
32 gross 2580 import matplotlib.pyplot as plt
33 lkettle 582
34 gross 2580
35 lkettle 582 ne=32 # number of cells in x_0 and x_1 directions
36     width=10000. # length in x_0 and x_1 directions
37 jgs 110 lam=3.462e9
38     mu=3.462e9
39     rho=1154.
40 gross 2580 tend=10. # to ran a full simulation change tend to 60.
41 gross 2513 alpha=0.7
42     t0=3.
43 jgs 108
44 gross 2580
45 gross 2513 U0=1. # maximum displacement
46     mkDir("data") # create directory data if it does not exist already.
47 jgs 110
48 gross 2513 def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
49 gross 2580 # lists to collect displacement at point source
50     ts, u_pc0,u_pc1,u_pc2=[], [], [], []
51 jgs 108 x=domain.getX()
52     # ... open new PDE ...
53 jgs 110 mypde=LinearPDE(domain)
54 gross 2474 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
55 lkettle 582 kronecker=identity(mypde.getDim())
56    
57 jfenwick 2455 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
58 lkettle 582
59 gross 2502 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
60 jgs 110 # ... set initial values ....
61 jgs 108 n=0
62 lkettle 582 # for first two time steps
63 gross 2502 u=Vector(0.,Solution(domain))
64     u_last=Vector(0.,Solution(domain))
65 jgs 110 t=0
66 lkettle 582
67     # define the location of the point source
68 gross 2502 L=Locator(domain,xc)
69 lkettle 582 # find potential at point source
70     u_pc=L.getValue(u)
71     print "u at point charge=",u_pc
72 gross 2580 ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
73 lkettle 582
74 jgs 108 while t<tend:
75 gross 2502 t+=h
76 jgs 108 # ... get current stress ....
77     g=grad(u)
78 lkettle 582 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
79 jgs 108 # ... get new acceleration ....
80 gross 2513 amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2*exp(1./2.-(t-t0)**2/alpha**2)
81 gross 2502 mypde.setValue(X=-stress, r=dunit*amplitude)
82 jgs 110 a=mypde.getSolution()
83     # ... get new displacement ...
84     u_new=2*u-u_last+h**2*a
85     # ... shift displacements ....
86     u_last=u
87     u=u_new
88 jgs 108 n+=1
89 jgs 110 print n,"-th time step t ",t
90 lkettle 582 u_pc=L.getValue(u)
91     print "u at point charge=",u_pc
92 gross 2580 ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
93 lkettle 582
94     # ... save current acceleration in units of gravity and displacements
95     if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
96 jongui 1161 displacement = length(u), tensor = stress, Ux = u[0] )
97 gross 2580 return ts, u_pc0,u_pc1,u_pc2
98 jgs 108
99 gross 2580 #
100     # create domain:
101     #
102 gross 2502 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
103 gross 2580 #
104     # sety time step size:
105     #
106 gross 2502 h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
107     print "time step size = ",h
108 gross 2580 #
109 gross 2513 # spherical source at middle of bottom face
110 gross 2580 #
111 gross 2513 xc=[width/2.,width/2.,0.]
112     # define small radius around point xc
113     src_radius = 0.03*width
114     print "src_radius = ",src_radius
115 gross 2580 #
116     # run it
117     #
118     ts, u_pc0,u_pc1,u_pc2 = wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)
119     #
120     # create a plot:
121     #
122     if getMPIRankWorld() == 0:
123     plt.title("Displacement at Point Source")
124     plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
125     plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
126     plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
127     plt.xlabel('time')
128     plt.ylabel('displacement')
129     plt.legend()
130     plt.savefig('u_pc.png', format='png')
131     # or save displacement
132     u_pc_data=FileWriter('./data/U_pc.out')
133     for i in xrange(len(ts)) :
134     u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
135     u_pc_data.close()
136 jgs 108

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/doc/examples/usersguide/wave.py:2682-2741 /branches/restext/doc/examples/usersguide/wave.py:2610-2624 /trunk/doc/examples/usersguide/wave.py:1388-2483

  ViewVC Help
Powered by ViewVC 1.1.26