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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2502 - (show annotations)
Tue Jun 30 05:49:22 2009 UTC (9 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 3576 byte(s)
tests for examples added. utest is not working yet.
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 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-2008 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.finley import Brick
26 from numpy import identity,zeros,ones
27
28 if not os.path.isdir("data"):
29 print "\nCreating subdirectory 'data'\n"
30 if getMPIRankWorld()==0: os.mkdir("data")
31
32 ne=32 # number of cells in x_0 and x_1 directions
33 width=10000. # length in x_0 and x_1 directions
34 lam=3.462e9
35 mu=3.462e9
36 rho=1154.
37 tend=0.05 # to ran a full simulation change tend to 60.
38 alpha=0.3
39
40 U0=0.01 # amplitude of point source
41
42 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
43 x=domain.getX()
44 # ... open new PDE ...
45 mypde=LinearPDE(domain)
46 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
47 kronecker=identity(mypde.getDim())
48
49 # spherical source at middle of bottom face
50
51 xc=[width/2.,width/2.,0.]
52 # define small radius around point xc
53 src_radius = 0.03*width
54 print "src_radius = ",src_radius
55
56 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
57
58 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
59 # ... set initial values ....
60 n=0
61 # initial value of displacement at point source is constant (U0=0.01)
62 # for first two time steps
63 u=Vector(0.,Solution(domain))
64 u_last=Vector(0.,Solution(domain))
65 t=0
66
67 # define the location of the point source
68 L=Locator(domain,xc)
69 # find potential at point source
70 u_pc=L.getValue(u)
71 print "u at point charge=",u_pc
72 # open file to save displacement at point source
73 u_pc_data=FileWriter('./data/U_pc.out')
74 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
75
76 while t<tend:
77 t+=h
78 # ... get current stress ....
79 g=grad(u)
80 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
81 # ... get new acceleration ....
82 amplitude=U0*2*exp(1)/alpha**2*(1-5*(t/alpha)**2+2*(t/alpha)**4)*exp(-(t/alpha)**2)
83 mypde.setValue(X=-stress, r=dunit*amplitude)
84 a=mypde.getSolution()
85 # ... get new displacement ...
86 u_new=2*u-u_last+h**2*a
87 # ... shift displacements ....
88 u_last=u
89 u=u_new
90 n+=1
91 print n,"-th time step t ",t
92 u_pc=L.getValue(u)
93 print "u at point charge=",u_pc
94 # save displacements at point source to file for t > 0
95 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
96
97 # ... save current acceleration in units of gravity and displacements
98 if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
99 displacement = length(u), tensor = stress, Ux = u[0] )
100
101 u_pc_data.close()
102
103 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
104 h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
105 print "time step size = ",h
106 wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
107

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26