| 1 |
# $Id$ |
# $Id$ |
| 2 |
from esys.escript import * |
from esys.escript import * |
| 3 |
import esys.finley |
from esys.escript.linearPDEs import LinearPDE |
| 4 |
import numarray |
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,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax): |
def wavePropagation(domain,h,tend,lam,mu,rho,s_tt): |
|
# ... get characteristic function of impact region: |
|
| 20 |
x=domain.getX() |
x=domain.getX() |
|
location=x[0].whereZero()*(length(x-xc)-r).whereNegative() |
|
| 21 |
# ... open new PDE ... |
# ... open new PDE ... |
| 22 |
myPDE=LinearPDE(mydomain) |
mypde=LinearPDE(domain) |
| 23 |
myPDE.setLumpingOn() |
mypde.setSolverMethod(mypde.LUMPING) |
| 24 |
myPDE.setValue(D=numarray.identity(myPDE.getDim())*rho,q=location*numarray.identity(myPDE.getDim())[:,0]) |
mypde.setValue(D=kronecker(mypde.getDim())*rho, \ |
| 25 |
# ... set the initial values : |
q=whereZero(x[0])*kronecker(mypde.getDim())[1,:]) |
| 26 |
t=dt |
# ... set initial values .... |
| 27 |
n=0 |
n=0 |
| 28 |
u=0 |
u=Vector(0,ContinuousFunction(domain)) |
| 29 |
v=0 |
u_last=Vector(0,ContinuousFunction(domain)) |
| 30 |
a=0 |
t=0 |
| 31 |
while t<tend: |
while t<tend: |
|
# ... up-date displacement .... |
|
|
u=u+dt*v+dt**2/2*a |
|
| 32 |
# ... get current stress .... |
# ... get current stress .... |
| 33 |
g=grad(u) |
g=grad(u) |
| 34 |
stress=lame_lambda*trace(g)+lame_mu*(g+transpose(g)) |
stress=lam*trace(g)*kronecker(mypde.getDim())+mu*(g+transpose(g)) |
| 35 |
# ... get new acceleration .... |
# ... get new acceleration .... |
| 36 |
myPDE.setValue(X=stress,q=impact_location, \ |
mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker(mypde.getDim())[1,:]) |
| 37 |
r=umax/tau**2*(6*t/tau-9*(t/tau)^4)*exp(-(t/tau)^3)*numarray.identity(myPDE.getDim())[:,0]) |
a=mypde.getSolution() |
| 38 |
a_new=myPDE.getSolution() |
# ... get new displacement ... |
| 39 |
# ... update velocity ... |
u_new=2*u-u_last+h**2*a |
| 40 |
v=v+h/2*(a+a_new) |
# ... shift displacements .... |
| 41 |
# ... next time step ... |
u_last=u |
| 42 |
a=a_new |
u=u_new |
| 43 |
t+=dt |
t+=h |
| 44 |
n+=1 |
n+=1 |
| 45 |
# ... save current displacement: |
print n,"-th time step t ",t |
| 46 |
if n%10: u.saveDX("u.%i.dx"%n) |
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 |
ne=6 |
mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width) |
| 52 |
lame_lambda=3.462e9 |
wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt) |
|
lame_mu=3.462e9 |
|
|
rho=1154. |
|
|
tau=2. |
|
|
umax=15. |
|
|
xc=[0,1000,1000] |
|
|
r=1. |
|
|
tend=10. |
|
|
dt=1./5.*sqrt(rho/(lame_lambda+2*lame_mu))(20000./ne) |
|
|
print "step size = ",dt |
|
|
mydomain=finely.Brick(ne,10*ne,10*ne,l0=20000,l1=200000,l2=200000) |
|
|
wavePropagation(mydomain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax) |
|
| 53 |
|
|
Parent Directory
| 
Revision Log
| 
Patch