/[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 6756 - (hide annotations)
Thu Nov 29 07:23:43 2018 UTC (7 months, 2 weeks ago) by aellery
File MIME type: text/x-python
File size: 5006 byte(s)
Temporarily undoing last commit.
1 jfenwick 3981 ##############################################################################
2 ksteube 1811 #
3 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 ksteube 1811 #
6     # Primary Business: Queensland, Australia
7 jfenwick 6112 # Licensed under the Apache License, version 2.0
8     # http://www.apache.org/licenses/LICENSE-2.0
9 ksteube 1811 #
10 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 sshaw 5707 from __future__ import division, print_function
16 ksteube 1811
17 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ksteube 1811 Primary Business: Queensland, Australia"""
20 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
21     http://www.apache.org/licenses/LICENSE-2.0"""
22 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
23 ksteube 1811
24 jfenwick 2600 import matplotlib
25     matplotlib.use('agg') #For interactive use, you can comment out this line
26     #It's just here to make testing easier
27 gross 2580 import matplotlib.pyplot as plt
28 caltinay 4088 from numpy import zeros,ones
29 jfenwick 6094 import numpy
30 caltinay 4087 from esys.escript import *
31 sshaw 4821 from esys.escript.linearPDEs import LinearPDE, SolverOptions
32 caltinay 4087 from esys.escript.pdetools import Locator
33 sshaw 5288 try:
34 sshaw 5292 from esys.dudley import Brick
35 sshaw 5288 HAVE_DUDLEY = True
36     except ImportError:
37     HAVE_DUDLEY = False
38 caltinay 4087 from esys.weipa import saveVTK
39 lkettle 582
40 sshaw 5288 if not HAVE_DUDLEY:
41     print("Dudley module not available")
42     else:
43     ne=32 # number of cells in x_0 and x_1 directions
44     width=10000. # length in x_0 and x_1 directions
45     lam=3.462e9
46     mu=3.462e9
47     rho=1154.
48     tend=10. # to ran a full simulation change tend to 60.
49     alpha=0.7
50     t0=3.
51 gross 2580
52 jgs 108
53 sshaw 5288 U0=1. # maximum displacement
54     mkDir("data") # create directory data if it does not exist already.
55 gross 2580
56 sshaw 5288 def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
57     # lists to collect displacement at point source
58     ts, u_pc0,u_pc1,u_pc2=[], [], [], []
59     x=domain.getX()
60     # ... open new PDE ...
61     mypde=LinearPDE(domain)
62     mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
63     kron=kronecker(mypde.getDim())
64 jgs 110
65 sshaw 5288 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
66 lkettle 582
67 sshaw 5288 mypde.setValue(D=kron*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
68     # ... set initial values ....
69     n=0
70     # for first two time steps
71     u=Vector(0.,Solution(domain))
72     u_last=Vector(0.,Solution(domain))
73     t=0
74 lkettle 582
75 sshaw 5288 # define the location of the point source
76     L=Locator(domain,xc)
77     # find potential at point source
78     u_pc=L.getValue(u)
79     print("u at point charge = %s"%u_pc)
80     ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
81    
82     while t<tend:
83     t+=h
84     # ... get current stress ....
85     g=grad(u)
86     stress=lam*trace(g)*kron+mu*(g+transpose(g))
87     # ... get new acceleration ....
88     amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2*exp(1./2.-(t-t0)**2/alpha**2)
89     mypde.setValue(X=-stress, r=dunit*amplitude)
90     a=mypde.getSolution()
91     # ... get new displacement ...
92     u_new=2*u-u_last+h**2*a
93     # ... shift displacements ....
94     u_last=u
95     u=u_new
96     n+=1
97     print("time step %d, t = %s"%(n,t))
98     u_pc=L.getValue(u)
99     print("u at point charge = %s"%u_pc)
100     ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
101    
102     # ... save current acceleration in units of gravity and displacements
103     if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
104     displacement = length(u), tensor = stress, Ux = u[0] )
105     return ts, u_pc0,u_pc1,u_pc2
106 lkettle 582
107 sshaw 5288 #
108     # create domain:
109     #
110     mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
111     #
112     # sety time step size:
113     #
114     h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
115     print("time step size = %s"%h)
116     #
117     # spherical source at middle of bottom face
118     #
119     xc=[width/2.,width/2.,0.]
120     # define small radius around point xc
121     src_radius = 0.03*width
122     print("src_radius = %s"%src_radius)
123     #
124     # run it
125     #
126     ts, u_pc0,u_pc1,u_pc2 = wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)
127     #
128     # create a plot:
129     #
130     if getMPIRankWorld() == 0:
131     plt.title("Displacement at Point Source")
132     plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
133     plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
134     plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
135     plt.xlabel('time')
136     plt.ylabel('displacement')
137     plt.legend()
138     plt.savefig('u_pc.png', format='png')
139     # or save displacement
140     u_pc_data=FileWriter('./data/U_pc.out')
141     for i in range(len(ts)) :
142     u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
143     u_pc_data.close()
144 jgs 108

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/doc/examples/usersguide/wave.py:5567-5588 /branches/lapack2681/doc/examples/usersguide/wave.py:2682-2741 /branches/pasowrap/doc/examples/usersguide/wave.py:3661-3674 /branches/py3_attempt2/doc/examples/usersguide/wave.py:3871-3891 /branches/restext/doc/examples/usersguide/wave.py:2610-2624 /branches/ripleygmg_from_3668/doc/examples/usersguide/wave.py:3669-3791 /branches/symbolic_from_3470/doc/examples/usersguide/wave.py:3471-3974 /branches/symbolic_from_3470/ripley/test/python/doc/examples/usersguide/wave.py:3517-3974 /branches/trilinos_from_5897/doc/examples/usersguide/wave.py:5898-6118 /release/4.0/doc/examples/usersguide/wave.py:5380-5406 /trunk/doc/examples/usersguide/wave.py:1388-2483 /trunk/ripley/test/python/doc/examples/usersguide/wave.py:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26