/[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 2549 - (hide annotations)
Mon Jul 20 06:43:47 2009 UTC (9 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 3476 byte(s)
Remainder of copyright date fixes
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 lkettle 582
28     ne=32 # number of cells in x_0 and x_1 directions
29     width=10000. # length in x_0 and x_1 directions
30 jgs 110 lam=3.462e9
31     mu=3.462e9
32     rho=1154.
33 gross 2513 tend=20. # to ran a full simulation change tend to 60.
34     alpha=0.7
35     t0=3.
36 jgs 108
37 gross 2513 U0=1. # maximum displacement
38     mkDir("data") # create directory data if it does not exist already.
39 jgs 110
40 gross 2513 def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
41 jgs 108 x=domain.getX()
42     # ... open new PDE ...
43 jgs 110 mypde=LinearPDE(domain)
44 gross 2474 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
45 lkettle 582 kronecker=identity(mypde.getDim())
46    
47 jfenwick 2455 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
48 lkettle 582
49 gross 2502 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
50 jgs 110 # ... set initial values ....
51 jgs 108 n=0
52 lkettle 582 # for first two time steps
53 gross 2502 u=Vector(0.,Solution(domain))
54     u_last=Vector(0.,Solution(domain))
55 jgs 110 t=0
56 lkettle 582
57     # define the location of the point source
58 gross 2502 L=Locator(domain,xc)
59 lkettle 582 # find potential at point source
60     u_pc=L.getValue(u)
61     print "u at point charge=",u_pc
62     # open file to save displacement at point source
63 gross 2420 u_pc_data=FileWriter('./data/U_pc.out')
64 gross 2502 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
65 lkettle 582
66 jgs 108 while t<tend:
67 gross 2502 t+=h
68 jgs 108 # ... get current stress ....
69     g=grad(u)
70 lkettle 582 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
71 jgs 108 # ... get new acceleration ....
72 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)
73 gross 2502 mypde.setValue(X=-stress, r=dunit*amplitude)
74 jgs 110 a=mypde.getSolution()
75     # ... get new displacement ...
76     u_new=2*u-u_last+h**2*a
77     # ... shift displacements ....
78     u_last=u
79     u=u_new
80 jgs 108 n+=1
81 jgs 110 print n,"-th time step t ",t
82 lkettle 582 u_pc=L.getValue(u)
83     print "u at point charge=",u_pc
84     # save displacements at point source to file for t > 0
85 gross 2502 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
86 lkettle 582
87     # ... save current acceleration in units of gravity and displacements
88     if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
89 jongui 1161 displacement = length(u), tensor = stress, Ux = u[0] )
90 jgs 108
91 lkettle 582 u_pc_data.close()
92    
93 gross 2502 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
94     h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
95     print "time step size = ",h
96 gross 2513 # spherical source at middle of bottom face
97     xc=[width/2.,width/2.,0.]
98     # define small radius around point xc
99     src_radius = 0.03*width
100     print "src_radius = ",src_radius
101     wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)
102 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