user/examples/wave.py

# $Id$
from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Brick
ne=5           # number of cells in x_0-direction
depth=10000.   # length in x_0-direction
width=100000.  # length in x_1 and x_2 direction
lam=3.462e9
mu=3.462e9
rho=1154.
tau=10.
umax=2.
tend=60
h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)
print "time step size = ",h

def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)

def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):
   x=domain.getX()
   # ... open new PDE ...
   mypde=LinearPDE(domain)
   mypde.setSolverMethod(mypde.LUMPING)
   mypde.setValue(D=kronecker(mypde.getDim())*rho, \
                  q=whereZero(x[0])*kronecker(mypde.getDim())[1,:])
   # ... set initial values ....
   n=0
   u=Vector(0,ContinuousFunction(domain))
   u_last=Vector(0,ContinuousFunction(domain))
   t=0
   while t<tend:
     # ... get current stress ....
     g=grad(u)
     stress=lam*trace(g)*kronecker(mypde.getDim())+mu*(g+transpose(g))
     # ... get new acceleration ....
     mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker(mypde.getDim())[1,:])
     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
     print "a=",inf(a),sup(a)
     print "u=",inf(u),sup(u)
     # ... save current acceleration in units of gravity
     if n%10==0: saveVTK("u.%i.xml"%(n/10),a=length(a)/9.81)

mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)
wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)

1	jgs	108	# $Id$
2			from esys.escript import *
3	gross	327	from esys.escript.linearPDEs import LinearPDE
4	jgs	110	from esys.finley import Brick
5	gross	327	ne=5 # number of cells in x_0-direction
6	jgs	110	depth=10000. # length in x_0-direction
7			width=100000. # length in x_1 and x_2 direction
8			lam=3.462e9
9			mu=3.462e9
10			rho=1154.
11			tau=10.
12			umax=2.
13			tend=60
14			h=1./5.sqrt(rho/(lam+2mu))*(depth/ne)
15			print "time step size = ",h
16	jgs	108
17	jgs	110	def s_tt(t): return umax/tau*2(6t/tau-9(t/tau)*4)exp(-(t/tau)**3)
18
19			def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):
20	jgs	108	x=domain.getX()
21			# ... open new PDE ...
22	jgs	110	mypde=LinearPDE(domain)
23	gross	327	mypde.setSolverMethod(mypde.LUMPING)
24			mypde.setValue(D=kronecker(mypde.getDim())*rho, \
25			q=whereZero(x[0])*kronecker(mypde.getDim())[1,:])
26	jgs	110	# ... set initial values ....
27	jgs	108	n=0
28	jgs	110	u=Vector(0,ContinuousFunction(domain))
29			u_last=Vector(0,ContinuousFunction(domain))
30			t=0
31	jgs	108	while t<tend:
32			# ... get current stress ....
33			g=grad(u)
34	gross	327	stress=lamtrace(g)kronecker(mypde.getDim())+mu*(g+transpose(g))
35	jgs	108	# ... get new acceleration ....
36	gross	327	mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker(mypde.getDim())[1,:])
37	jgs	110	a=mypde.getSolution()
38			# ... get new displacement ...
39			u_new=2u-u_last+h2a
40			# ... shift displacements ....
41			u_last=u
42			u=u_new
43			t+=h
44	jgs	108	n+=1
45	jgs	110	print n,"-th time step t ",t
46			print "a=",inf(a),sup(a)
47			print "u=",inf(u),sup(u)
48			# ... save current acceleration in units of gravity
49	gross	327	if n%10==0: saveVTK("u.%i.xml"%(n/10),a=length(a)/9.81)
50	jgs	108
51	jgs	110	mydomain=Brick(ne,int(width/depth)ne,int(width/depth)ne,l0=depth,l1=width,l2=width)
52			wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
53	jgs	108
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision