/[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 2502 - (hide annotations)
Tue Jun 30 05:49:22 2009 UTC (10 years, 4 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 ksteube 1215
2 ksteube 1811 ########################################################
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 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 lkettle 582
28 ksteube 1215 if not os.path.isdir("data"):
29     print "\nCreating subdirectory 'data'\n"
30 gross 2502 if getMPIRankWorld()==0: os.mkdir("data")
31 ksteube 1215
32 lkettle 582 ne=32 # number of cells in x_0 and x_1 directions
33     width=10000. # length in x_0 and x_1 directions
34 jgs 110 lam=3.462e9
35     mu=3.462e9
36     rho=1154.
37 gross 2502 tend=0.05 # to ran a full simulation change tend to 60.
38     alpha=0.3
39 jgs 108
40 lkettle 582 U0=0.01 # amplitude of point source
41 jgs 110
42 lkettle 582 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
43 jgs 108 x=domain.getX()
44     # ... open new PDE ...
45 jgs 110 mypde=LinearPDE(domain)
46 gross 2474 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
47 lkettle 582 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 gross 2502 src_radius = 0.03*width
54 lkettle 582 print "src_radius = ",src_radius
55    
56 jfenwick 2455 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
57 lkettle 582
58 gross 2502 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
59 jgs 110 # ... set initial values ....
60 jgs 108 n=0
61 lkettle 582 # initial value of displacement at point source is constant (U0=0.01)
62     # 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     # open file to save displacement at point source
73 gross 2420 u_pc_data=FileWriter('./data/U_pc.out')
74 gross 2502 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
75 lkettle 582
76 jgs 108 while t<tend:
77 gross 2502 t+=h
78 jgs 108 # ... get current stress ....
79     g=grad(u)
80 lkettle 582 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
81 jgs 108 # ... get new acceleration ....
82 gross 2502 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 jgs 110 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 jgs 108 n+=1
91 jgs 110 print n,"-th time step t ",t
92 lkettle 582 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 gross 2502 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
96 lkettle 582
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 jongui 1161 displacement = length(u), tensor = stress, Ux = u[0] )
100 jgs 108
101 lkettle 582 u_pc_data.close()
102    
103 gross 2502 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 lkettle 582 wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
107 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