/[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 3675 - (hide annotations)
Thu Nov 17 00:53:38 2011 UTC (7 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 4309 byte(s)
pasowrap joins the trunk.

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