1 |
# $Id$ |
2 |
from esys.escript import * |
3 |
from esys.escript.linearPDEs import LinearPDE |
4 |
from esys.finley import Brick |
5 |
ne=5 # number of cells in x_0-direction |
6 |
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+2*mu))*(depth/ne) |
15 |
print "time step size = ",h |
16 |
|
17 |
def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3) |
18 |
|
19 |
def wavePropagation(domain,h,tend,lam,mu,rho,s_tt): |
20 |
x=domain.getX() |
21 |
# ... open new PDE ... |
22 |
mypde=LinearPDE(domain) |
23 |
mypde.setSolverMethod(mypde.LUMPING) |
24 |
mypde.setValue(D=kronecker(mypde.getDim())*rho, \ |
25 |
q=whereZero(x[0])*kronecker(mypde.getDim())[1,:]) |
26 |
# ... set initial values .... |
27 |
n=0 |
28 |
u=Vector(0,ContinuousFunction(domain)) |
29 |
u_last=Vector(0,ContinuousFunction(domain)) |
30 |
t=0 |
31 |
while t<tend: |
32 |
# ... get current stress .... |
33 |
g=grad(u) |
34 |
stress=lam*trace(g)*kronecker(mypde.getDim())+mu*(g+transpose(g)) |
35 |
# ... get new acceleration .... |
36 |
mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker(mypde.getDim())[1,:]) |
37 |
a=mypde.getSolution() |
38 |
# ... get new displacement ... |
39 |
u_new=2*u-u_last+h**2*a |
40 |
# ... shift displacements .... |
41 |
u_last=u |
42 |
u=u_new |
43 |
t+=h |
44 |
n+=1 |
45 |
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 |
if n%10==0: saveVTK("u.%i.xml"%(n/10),a=length(a)/9.81) |
50 |
|
51 |
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 |
|