/[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 2548 - (show annotations)
Mon Jul 20 06:20:06 2009 UTC (9 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 3476 byte(s)
Updating copyright notices
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 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 ne=32 # number of cells in x_0 and x_1 directions
29 width=10000. # length in x_0 and x_1 directions
30 lam=3.462e9
31 mu=3.462e9
32 rho=1154.
33 tend=20. # to ran a full simulation change tend to 60.
34 alpha=0.7
35 t0=3.
36
37 U0=1. # maximum displacement
38 mkDir("data") # create directory data if it does not exist already.
39
40 def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
41 x=domain.getX()
42 # ... open new PDE ...
43 mypde=LinearPDE(domain)
44 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
45 kronecker=identity(mypde.getDim())
46
47 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
48
49 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
50 # ... set initial values ....
51 n=0
52 # for first two time steps
53 u=Vector(0.,Solution(domain))
54 u_last=Vector(0.,Solution(domain))
55 t=0
56
57 # define the location of the point source
58 L=Locator(domain,xc)
59 # 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 u_pc_data=FileWriter('./data/U_pc.out')
64 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
65
66 while t<tend:
67 t+=h
68 # ... get current stress ....
69 g=grad(u)
70 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
71 # ... get new acceleration ....
72 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 mypde.setValue(X=-stress, r=dunit*amplitude)
74 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 n+=1
81 print n,"-th time step t ",t
82 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 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
86
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 displacement = length(u), tensor = stress, Ux = u[0] )
90
91 u_pc_data.close()
92
93 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 # 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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26