| 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 |
|
Parent Directory
| 
Revision Log