/[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 2580 - (hide annotations)
Tue Aug 4 08:24:12 2009 UTC (9 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 4138 byte(s)
another example for the usage of matplot lib added to users guide.
1 ksteube 1215
2 ksteube 1811 ########################################################
3     #
4 jfenwick 2548 # Copyright (c) 2003-2009 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 2549 __copyright__="""Copyright (c) 2003-2009 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 gross 2580 import matplotlib.pyplot as plt
28 lkettle 582
29 gross 2580
30 lkettle 582 ne=32 # number of cells in x_0 and x_1 directions
31     width=10000. # length in x_0 and x_1 directions
32 jgs 110 lam=3.462e9
33     mu=3.462e9
34     rho=1154.
35 gross 2580 tend=10. # to ran a full simulation change tend to 60.
36 gross 2513 alpha=0.7
37     t0=3.
38 jgs 108
39 gross 2580
40 gross 2513 U0=1. # maximum displacement
41     mkDir("data") # create directory data if it does not exist already.
42 jgs 110
43 gross 2513 def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
44 gross 2580 # lists to collect displacement at point source
45     ts, u_pc0,u_pc1,u_pc2=[], [], [], []
46 jgs 108 x=domain.getX()
47     # ... open new PDE ...
48 jgs 110 mypde=LinearPDE(domain)
49 gross 2474 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
50 lkettle 582 kronecker=identity(mypde.getDim())
51    
52 jfenwick 2455 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
53 lkettle 582
54 gross 2502 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
55 jgs 110 # ... set initial values ....
56 jgs 108 n=0
57 lkettle 582 # for first two time steps
58 gross 2502 u=Vector(0.,Solution(domain))
59     u_last=Vector(0.,Solution(domain))
60 jgs 110 t=0
61 lkettle 582
62     # define the location of the point source
63 gross 2502 L=Locator(domain,xc)
64 lkettle 582 # find potential at point source
65     u_pc=L.getValue(u)
66     print "u at point charge=",u_pc
67 gross 2580 ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
68 lkettle 582
69 jgs 108 while t<tend:
70 gross 2502 t+=h
71 jgs 108 # ... get current stress ....
72     g=grad(u)
73 lkettle 582 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
74 jgs 108 # ... get new acceleration ....
75 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)
76 gross 2502 mypde.setValue(X=-stress, r=dunit*amplitude)
77 jgs 110 a=mypde.getSolution()
78     # ... get new displacement ...
79     u_new=2*u-u_last+h**2*a
80     # ... shift displacements ....
81     u_last=u
82     u=u_new
83 jgs 108 n+=1
84 jgs 110 print n,"-th time step t ",t
85 lkettle 582 u_pc=L.getValue(u)
86     print "u at point charge=",u_pc
87 gross 2580 ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
88 lkettle 582
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 jongui 1161 displacement = length(u), tensor = stress, Ux = u[0] )
92 gross 2580 return ts, u_pc0,u_pc1,u_pc2
93 jgs 108
94 gross 2580 #
95     # create domain:
96     #
97 gross 2502 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
98 gross 2580 #
99     # sety time step size:
100     #
101 gross 2502 h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
102     print "time step size = ",h
103 gross 2580 #
104 gross 2513 # spherical source at middle of bottom face
105 gross 2580 #
106 gross 2513 xc=[width/2.,width/2.,0.]
107     # define small radius around point xc
108     src_radius = 0.03*width
109     print "src_radius = ",src_radius
110 gross 2580 #
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 jgs 108

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo

  ViewVC Help
Powered by ViewVC 1.1.26