examples/usersguide/wave.py


########################################################
#
# Copyright (c) 2003-2010 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2010 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

from esys.escript import *
from esys.escript.pdetools import Locator
from esys.escript.linearPDEs import LinearPDE
from esys.dudley import Brick
from esys.weipa import saveVTK
from numpy import identity,zeros,ones
import matplotlib

matplotlib.use('agg')    #For interactive use, you can comment out this line
#It's just here to make testing easier

import matplotlib.pyplot as plt


ne=32          # number of cells in x_0 and x_1 directions
width=10000.  # length in x_0 and x_1 directions
lam=3.462e9
mu=3.462e9
rho=1154.
tend=10. # to ran a full simulation change tend to 60.
alpha=0.7
t0=3.


U0=1. # maximum displacement
mkDir("data") # create directory data if it does not exist already.

def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
   # lists to collect displacement at point source
   ts, u_pc0,u_pc1,u_pc2=[], [], [], []
   x=domain.getX()
   # ... open new PDE ...
   mypde=LinearPDE(domain)
   mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING)
   kronecker=identity(mypde.getDim())

   dunit=numpy.array([1.,0.,0.]) # defines direction of point source

   mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
   # ... set initial values ....
   n=0
   # for first two time steps
   u=Vector(0.,Solution(domain))
   u_last=Vector(0.,Solution(domain))
   t=0

   # define the location of the point source 
   L=Locator(domain,xc)
   # find potential at point source
   u_pc=L.getValue(u)
   print "u at point charge=",u_pc
   ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
 
   while t<tend:
     t+=h
     # ... get current stress ....
     g=grad(u)
     stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
     # ... get new acceleration ....
     amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2*exp(1./2.-(t-t0)**2/alpha**2)
     mypde.setValue(X=-stress, r=dunit*amplitude)
     a=mypde.getSolution()
     # ... get new displacement ...
     u_new=2*u-u_last+h**2*a
     # ... shift displacements ....
     u_last=u
     u=u_new
     n+=1
     print n,"-th time step t ",t
     u_pc=L.getValue(u)
     print "u at point charge=",u_pc
     ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
 
     # ... save current acceleration in units of gravity and displacements 
     if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
     displacement = length(u), tensor = stress, Ux = u[0] )
   return ts, u_pc0,u_pc1,u_pc2

#
# create domain:
#
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
#
#  sety time step size:
#
h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
print "time step size = ",h
#
#  spherical source at middle of bottom face
#
xc=[width/2.,width/2.,0.]
# define small radius around point xc
src_radius = 0.03*width
print "src_radius = ",src_radius
#
# run it
#
ts, u_pc0,u_pc1,u_pc2 = wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)
#
# create a plot:
#
if getMPIRankWorld() == 0:
    plt.title("Displacement at Point Source")
    plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
    plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
    plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
    plt.xlabel('time')
    plt.ylabel('displacement')
    plt.legend()
    plt.savefig('u_pc.png', format='png')
# or save displacement
u_pc_data=FileWriter('./data/U_pc.out')
for i in xrange(len(ts)) :
    u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
u_pc_data.close()

1	ksteube	1215
2	ksteube	1811	########################################################
3			#
4	jfenwick	2881	# Copyright (c) 2003-2010 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	2881	__copyright__="""Copyright (c) 2003-2010 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	jfenwick	3259	from esys.dudley import Brick
26	caltinay	3346	from esys.weipa import saveVTK
27	jfenwick	2455	from numpy import identity,zeros,ones
28	jfenwick	2600	import matplotlib
29
30			matplotlib.use('agg') #For interactive use, you can comment out this line
31			#It's just here to make testing easier
32
33	gross	2580	import matplotlib.pyplot as plt
34	lkettle	582
35	gross	2580
36	lkettle	582	ne=32 # number of cells in x_0 and x_1 directions
37			width=10000. # length in x_0 and x_1 directions
38	jgs	110	lam=3.462e9
39			mu=3.462e9
40			rho=1154.
41	gross	2580	tend=10. # to ran a full simulation change tend to 60.
42	gross	2513	alpha=0.7
43			t0=3.
44	jgs	108
45	gross	2580
46	gross	2513	U0=1. # maximum displacement
47			mkDir("data") # create directory data if it does not exist already.
48	jgs	110
49	gross	2513	def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
50	gross	2580	# lists to collect displacement at point source
51			ts, u_pc0,u_pc1,u_pc2=[], [], [], []
52	jgs	108	x=domain.getX()
53			# ... open new PDE ...
54	jgs	110	mypde=LinearPDE(domain)
55	gross	3379	mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING)
56	lkettle	582	kronecker=identity(mypde.getDim())
57
58	jfenwick	2455	dunit=numpy.array([1.,0.,0.]) # defines direction of point source
59	lkettle	582
60	gross	2502	mypde.setValue(D=kroneckerrho, q=whereNegative(length(x-xc)-src_radius)dunit)
61	jgs	110	# ... set initial values ....
62	jgs	108	n=0
63	lkettle	582	# for first two time steps
64	gross	2502	u=Vector(0.,Solution(domain))
65			u_last=Vector(0.,Solution(domain))
66	jgs	110	t=0
67	lkettle	582
68			# define the location of the point source
69	gross	2502	L=Locator(domain,xc)
70	lkettle	582	# find potential at point source
71			u_pc=L.getValue(u)
72			print "u at point charge=",u_pc
73	gross	2580	ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
74	lkettle	582
75	jgs	108	while t<tend:
76	gross	2502	t+=h
77	jgs	108	# ... get current stress ....
78			g=grad(u)
79	lkettle	582	stress=lamtrace(g)kronecker+mu*(g+transpose(g))
80	jgs	108	# ... get new acceleration ....
81	gross	2513	amplitude=U0(4(t-t0)3/alpha3-6(t-t0)/alpha)sqrt(2.)/alpha*2exp(1./2.-(t-t0)2/alpha2)
82	gross	2502	mypde.setValue(X=-stress, r=dunit*amplitude)
83	jgs	110	a=mypde.getSolution()
84			# ... get new displacement ...
85			u_new=2u-u_last+h2a
86			# ... shift displacements ....
87			u_last=u
88			u=u_new
89	jgs	108	n+=1
90	jgs	110	print n,"-th time step t ",t
91	lkettle	582	u_pc=L.getValue(u)
92			print "u at point charge=",u_pc
93	gross	2580	ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
94	lkettle	582
95			# ... save current acceleration in units of gravity and displacements
96			if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
97	jongui	1161	displacement = length(u), tensor = stress, Ux = u[0] )
98	gross	2580	return ts, u_pc0,u_pc1,u_pc2
99	jgs	108
100	gross	2580	#
101			# create domain:
102			#
103	gross	2502	mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
104	gross	2580	#
105			# sety time step size:
106			#
107	gross	2502	h=inf(1./5.)inf(sqrt(rho/(lam+2mu))*mydomain.getSize())
108			print "time step size = ",h
109	gross	2580	#
110	gross	2513	# spherical source at middle of bottom face
111	gross	2580	#
112	gross	2513	xc=[width/2.,width/2.,0.]
113			# define small radius around point xc
114			src_radius = 0.03*width
115			print "src_radius = ",src_radius
116	gross	2580	#
117			# run it
118			#
119			ts, u_pc0,u_pc1,u_pc2 = wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)
120			#
121			# create a plot:
122			#
123			if getMPIRankWorld() == 0:
124			plt.title("Displacement at Point Source")
125			plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
126			plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
127			plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
128			plt.xlabel('time')
129			plt.ylabel('displacement')
130			plt.legend()
131			plt.savefig('u_pc.png', format='png')
132			# or save displacement
133			u_pc_data=FileWriter('./data/U_pc.out')
134			for i in xrange(len(ts)) :
135			u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
136			u_pc_data.close()
137	jgs	108
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision
svn:mergeinfo	/branches/lapack2681/doc/examples/usersguide/wave.py:2682-2741 /branches/pasowrap/doc/examples/usersguide/wave.py:3661-3674 /branches/restext/doc/examples/usersguide/wave.py:2610-2624 /trunk/doc/examples/usersguide/wave.py:1388-2483