doc/examples/wave.py

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

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), 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	jgs	108	from esys.escript import *
2	lkettle	582	from esys.escript.pdetools import Locator
3	gross	327	from esys.escript.linearPDEs import LinearPDE
4	jgs	110	from esys.finley import Brick
5	lkettle	582	from numarray import identity,zeros,ones
6
7			ne=32 # number of cells in x_0 and x_1 directions
8			width=10000. # length in x_0 and x_1 directions
9	jgs	110	lam=3.462e9
10			mu=3.462e9
11			rho=1154.
12			tend=60
13	lkettle	582	h=(1./5.)sqrt(rho/(lam+2mu))*(width/ne)
14	jgs	110	print "time step size = ",h
15	jgs	108
16	lkettle	582	U0=0.01 # amplitude of point source
17	jgs	110
18	lkettle	582	def wavePropagation(domain,h,tend,lam,mu,rho,U0):
19	jgs	108	x=domain.getX()
20			# ... open new PDE ...
21	jgs	110	mypde=LinearPDE(domain)
22	lkettle	582	mypde.setSolverMethod(LinearPDE.LUMPING)
23			kronecker=identity(mypde.getDim())
24
25			# spherical source at middle of bottom face
26
27			xc=[width/2.,width/2.,0.]
28			# define small radius around point xc
29			# Lsup(x) returns the maximum value of the argument x
30			src_radius = 0.1*Lsup(domain.getSize())
31			print "src_radius = ",src_radius
32
33			dunit=numarray.array([1.,0.,0.]) # defines direction of point source
34
35			mypde.setValue(D=kronecker*rho)
36	jgs	110	# ... set initial values ....
37	jgs	108	n=0
38	lkettle	582	# initial value of displacement at point source is constant (U0=0.01)
39			# for first two time steps
40			u=U0whereNegative(length(x-xc)-src_radius)dunit
41			u_last=U0whereNegative(length(x-xc)-src_radius)dunit
42	jgs	110	t=0
43	lkettle	582
44			# define the location of the point source
45	jongui	1058	L=Locator(domain,numarray.array(xc))
46	lkettle	582	# find potential at point source
47			u_pc=L.getValue(u)
48			print "u at point charge=",u_pc
49
50			u_pc_x = u_pc[0]
51			u_pc_y = u_pc[1]
52			u_pc_z = u_pc[2]
53
54			# open file to save displacement at point source
55			u_pc_data=open('./data/U_pc.out','w')
56			u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
57
58	jgs	108	while t<tend:
59			# ... get current stress ....
60			g=grad(u)
61	lkettle	582	stress=lamtrace(g)kronecker+mu*(g+transpose(g))
62	jgs	108	# ... get new acceleration ....
63	lkettle	582	mypde.setValue(X=-stress)
64	jgs	110	a=mypde.getSolution()
65			# ... get new displacement ...
66			u_new=2u-u_last+h2a
67			# ... shift displacements ....
68			u_last=u
69			u=u_new
70			t+=h
71	jgs	108	n+=1
72	jgs	110	print n,"-th time step t ",t
73	lkettle	582	u_pc=L.getValue(u)
74			print "u at point charge=",u_pc
75
76			u_pc_x=u_pc[0]
77			u_pc_y=u_pc[1]
78			u_pc_z=u_pc[2]
79
80			# save displacements at point source to file for t > 0
81			u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
82
83			# ... save current acceleration in units of gravity and displacements
84			if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
85			displacement = length(u), Ux = u[0] )
86	jgs	108
87	lkettle	582	u_pc_data.close()
88
89			mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
90			wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
91	jgs	108
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision