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

Contents of /branches/symbolic_from_3470/doc/examples/usersguide/wave.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 4309 byte(s)
Fast forward to latest trunk revision 3788.

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 by University of Queensland
5 # 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 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 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 __url__="https://launchpad.net/escript-finley"
21
22 from esys.escript import *
23 from esys.escript.pdetools import Locator
24 from esys.escript.linearPDEs import LinearPDE
25 from esys.dudley import Brick
26 from esys.weipa import saveVTK
27 from numpy import identity,zeros,ones
28 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 import matplotlib.pyplot as plt
34
35
36 ne=32 # number of cells in x_0 and x_1 directions
37 width=10000. # length in x_0 and x_1 directions
38 lam=3.462e9
39 mu=3.462e9
40 rho=1154.
41 tend=10. # to ran a full simulation change tend to 60.
42 alpha=0.7
43 t0=3.
44
45
46 U0=1. # maximum displacement
47 mkDir("data") # create directory data if it does not exist already.
48
49 def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
50 # lists to collect displacement at point source
51 ts, u_pc0,u_pc1,u_pc2=[], [], [], []
52 x=domain.getX()
53 # ... open new PDE ...
54 mypde=LinearPDE(domain)
55 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING)
56 kronecker=identity(mypde.getDim())
57
58 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
59
60 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
61 # ... set initial values ....
62 n=0
63 # for first two time steps
64 u=Vector(0.,Solution(domain))
65 u_last=Vector(0.,Solution(domain))
66 t=0
67
68 # define the location of the point source
69 L=Locator(domain,xc)
70 # find potential at point source
71 u_pc=L.getValue(u)
72 print "u at point charge=",u_pc
73 ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
74
75 while t<tend:
76 t+=h
77 # ... get current stress ....
78 g=grad(u)
79 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
80 # ... get new acceleration ....
81 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 mypde.setValue(X=-stress, r=dunit*amplitude)
83 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 n+=1
90 print n,"-th time step t ",t
91 u_pc=L.getValue(u)
92 print "u at point charge=",u_pc
93 ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
94
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 displacement = length(u), tensor = stress, Ux = u[0] )
98 return ts, u_pc0,u_pc1,u_pc2
99
100 #
101 # create domain:
102 #
103 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
104 #
105 # sety time step size:
106 #
107 h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
108 print "time step size = ",h
109 #
110 # spherical source at middle of bottom face
111 #
112 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 #
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

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,3471-3788

  ViewVC Help
Powered by ViewVC 1.1.26