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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2420 - (show annotations)
Thu May 14 02:28:58 2009 UTC (10 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 3638 byte(s)
vtk writer supports now meta data and meta data schema
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 # You can shorten the execution time by reducing variable tend from 60 to 0.5
23
24 from esys.escript import *
25 from esys.escript.pdetools import Locator
26 from esys.escript.linearPDEs import LinearPDE
27 from esys.finley import Brick
28 from numarray import identity,zeros,ones
29
30 if not os.path.isdir("data"):
31 print "\nCreating subdirectory 'data'\n"
32 os.mkdir("data")
33
34 ne=32 # number of cells in x_0 and x_1 directions
35 width=10000. # length in x_0 and x_1 directions
36 lam=3.462e9
37 mu=3.462e9
38 rho=1154.
39 tend=60
40 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
41 print "time step size = ",h
42
43 U0=0.01 # amplitude of point source
44
45 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
46 x=domain.getX()
47 # ... open new PDE ...
48 mypde=LinearPDE(domain)
49 mypde.setSolverMethod(LinearPDE.LUMPING)
50 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 dunit=numarray.array([1.,0.,0.]) # defines direction of point source
61
62 mypde.setValue(D=kronecker*rho)
63 # ... set initial values ....
64 n=0
65 # 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 t=0
70
71 # define the location of the point source
72 L=Locator(domain,numarray.array(xc))
73 # 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 u_pc_data=FileWriter('./data/U_pc.out')
83 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
84
85 while t<tend:
86 # ... get current stress ....
87 g=grad(u)
88 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
89 # ... get new acceleration ....
90 mypde.setValue(X=-stress)
91 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 n+=1
99 print n,"-th time step t ",t
100 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 displacement = length(u), tensor = stress, Ux = u[0] )
113
114 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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26