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