doc/examples/wave.py


########################################################
#
# Copyright (c) 2003-2008 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-2008 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__="http://www.uq.edu.au/esscc/escript-finley"

# You can shorten the execution time by reducing variable tend from 60 to 0.5

from esys.escript import *
from esys.escript.pdetools import Locator
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Brick
from numarray import identity,zeros,ones

if not os.path.isdir("data"):
   print "\nCreating subdirectory 'data'\n"
   os.mkdir("data")

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=60
h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
print "time step size = ",h

U0=0.01 # amplitude of point source

def wavePropagation(domain,h,tend,lam,mu,rho,U0):
   x=domain.getX()
   # ... open new PDE ...
   mypde=LinearPDE(domain)
   mypde.setSolverMethod(LinearPDE.LUMPING)
   kronecker=identity(mypde.getDim())

   #  spherical source at middle of bottom face

   xc=[width/2.,width/2.,0.]
   # define small radius around point xc
   # Lsup(x) returns the maximum value of the argument x
   src_radius = 0.1*Lsup(domain.getSize())
   print "src_radius = ",src_radius

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

   mypde.setValue(D=kronecker*rho)
   # ... set initial values ....
   n=0
   # initial value of displacement at point source is constant (U0=0.01)
   # for first two time steps
   u=U0*whereNegative(length(x-xc)-src_radius)*dunit
   u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
   t=0

   # define the location of the point source 
   L=Locator(domain,numarray.array(xc))
   # find potential at point source
   u_pc=L.getValue(u)
   print "u at point charge=",u_pc
  
   u_pc_x = u_pc[0]
   u_pc_y = u_pc[1]
   u_pc_z = u_pc[2]

   # open file to save displacement at point source
   u_pc_data=open('./data/U_pc.out','w')
   u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
 
   while t<tend:
     # ... get current stress ....
     g=grad(u)
     stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
     # ... get new acceleration ....
     mypde.setValue(X=-stress)          
     a=mypde.getSolution()
     # ... get new displacement ...
     u_new=2*u-u_last+h**2*a
     # ... shift displacements ....
     u_last=u
     u=u_new
     t+=h
     n+=1
     print n,"-th time step t ",t
     u_pc=L.getValue(u)
     print "u at point charge=",u_pc
     
     u_pc_x=u_pc[0]
     u_pc_y=u_pc[1]
     u_pc_z=u_pc[2]
      
     # save displacements at point source to file for t > 0
     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
 
     # ... 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] )

   u_pc_data.close()
  
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
wavePropagation(mydomain,h,tend,lam,mu,rho,U0)

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__="http://www.uq.edu.au/esscc/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+2mu))*(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=U0whereNegative(length(x-xc)-src_radius)dunit
68	u_last=U0whereNegative(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=open('./data/U_pc.out','w')
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=lamtrace(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=2u-u_last+h2a
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
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision