/[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 2474 - (hide annotations)
Tue Jun 16 06:32:15 2009 UTC (9 years, 8 months ago) by gross
Original Path: trunk/doc/examples/wave.py
File MIME type: text/x-python
File size: 3663 byte(s)
linearPDEs has is now using the SolverOptions class to talk to PASO
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 ksteube 1215 # You can shorten the execution time by reducing variable tend from 60 to 0.5
23    
24 jgs 108 from esys.escript import *
25 lkettle 582 from esys.escript.pdetools import Locator
26 gross 327 from esys.escript.linearPDEs import LinearPDE
27 jgs 110 from esys.finley import Brick
28 jfenwick 2455 from numpy import identity,zeros,ones
29 lkettle 582
30 ksteube 1215 if not os.path.isdir("data"):
31     print "\nCreating subdirectory 'data'\n"
32     os.mkdir("data")
33    
34 lkettle 582 ne=32 # number of cells in x_0 and x_1 directions
35     width=10000. # length in x_0 and x_1 directions
36 jgs 110 lam=3.462e9
37     mu=3.462e9
38     rho=1154.
39     tend=60
40 lkettle 582 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
41 jgs 110 print "time step size = ",h
42 jgs 108
43 lkettle 582 U0=0.01 # amplitude of point source
44 jgs 110
45 lkettle 582 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
46 jgs 108 x=domain.getX()
47     # ... open new PDE ...
48 jgs 110 mypde=LinearPDE(domain)
49 gross 2474 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
50 lkettle 582 kronecker=identity(mypde.getDim())
51    
52     # spherical source at middle of bottom face
53    
54     xc=[width/2.,width/2.,0.]
55     # define small radius around point xc
56     # Lsup(x) returns the maximum value of the argument x
57     src_radius = 0.1*Lsup(domain.getSize())
58     print "src_radius = ",src_radius
59    
60 jfenwick 2455 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
61 lkettle 582
62     mypde.setValue(D=kronecker*rho)
63 jgs 110 # ... set initial values ....
64 jgs 108 n=0
65 lkettle 582 # initial value of displacement at point source is constant (U0=0.01)
66     # for first two time steps
67     u=U0*whereNegative(length(x-xc)-src_radius)*dunit
68     u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
69 jgs 110 t=0
70 lkettle 582
71     # define the location of the point source
72 jfenwick 2455 L=Locator(domain,numpy.array(xc))
73 lkettle 582 # find potential at point source
74     u_pc=L.getValue(u)
75     print "u at point charge=",u_pc
76    
77     u_pc_x = u_pc[0]
78     u_pc_y = u_pc[1]
79     u_pc_z = u_pc[2]
80    
81     # open file to save displacement at point source
82 gross 2420 u_pc_data=FileWriter('./data/U_pc.out')
83 lkettle 582 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
84    
85 jgs 108 while t<tend:
86     # ... get current stress ....
87     g=grad(u)
88 lkettle 582 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
89 jgs 108 # ... get new acceleration ....
90 lkettle 582 mypde.setValue(X=-stress)
91 jgs 110 a=mypde.getSolution()
92     # ... get new displacement ...
93     u_new=2*u-u_last+h**2*a
94     # ... shift displacements ....
95     u_last=u
96     u=u_new
97     t+=h
98 jgs 108 n+=1
99 jgs 110 print n,"-th time step t ",t
100 lkettle 582 u_pc=L.getValue(u)
101     print "u at point charge=",u_pc
102    
103     u_pc_x=u_pc[0]
104     u_pc_y=u_pc[1]
105     u_pc_z=u_pc[2]
106    
107     # save displacements at point source to file for t > 0
108     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
109    
110     # ... save current acceleration in units of gravity and displacements
111     if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
112 jongui 1161 displacement = length(u), tensor = stress, Ux = u[0] )
113 jgs 108
114 lkettle 582 u_pc_data.close()
115    
116     mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
117     wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
118 jgs 108

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26