# Diff of /trunk/doc/examples/usersguide/wave.py

revision 2484 by gross, Mon Jun 22 04:22:19 2009 UTC revision 2502 by gross, Tue Jun 30 05:49:22 2009 UTC
21
# You can shorten the execution time by reducing variable tend from 60 to 0.5

22  from esys.escript import *  from esys.escript import *
23  from esys.escript.pdetools import Locator  from esys.escript.pdetools import Locator
24  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
# Line 29  from numpy import identity,zeros,ones Line 27  from numpy import identity,zeros,ones
27
28  if not os.path.isdir("data"):  if not os.path.isdir("data"):
29     print "\nCreating subdirectory 'data'\n"     print "\nCreating subdirectory 'data'\n"
30     os.mkdir("data")     if getMPIRankWorld()==0:  os.mkdir("data")
31
32  ne=32          # number of cells in x_0 and x_1 directions  ne=32          # number of cells in x_0 and x_1 directions
33  width=10000.  # length in x_0 and x_1 directions  width=10000.  # length in x_0 and x_1 directions
34  lam=3.462e9  lam=3.462e9
35  mu=3.462e9  mu=3.462e9
36  rho=1154.  rho=1154.
37  tend=60  tend=0.05     # to ran a full simulation change tend to 60.
38  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)  alpha=0.3
print "time step size = ",h
39
40  U0=0.01 # amplitude of point source  U0=0.01 # amplitude of point source
41
# Line 53  def wavePropagation(domain,h,tend,lam,mu Line 50  def wavePropagation(domain,h,tend,lam,mu
50
51     xc=[width/2.,width/2.,0.]     xc=[width/2.,width/2.,0.]
52     # define small radius around point xc     # define small radius around point xc
53     # Lsup(x) returns the maximum value of the argument x     src_radius = 0.03*width
55
56     dunit=numpy.array([1.,0.,0.]) # defines direction of point source     dunit=numpy.array([1.,0.,0.]) # defines direction of point source
57
59     # ... set initial values ....     # ... set initial values ....
60     n=0     n=0
61     # initial value of displacement at point source is constant (U0=0.01)     # initial value of displacement at point source is constant (U0=0.01)
62     # for first two time steps     # for first two time steps
65     t=0     t=0
66
67     # define the location of the point source     # define the location of the point source
68     L=Locator(domain,numpy.array(xc))     L=Locator(domain,xc)
69     # find potential at point source     # find potential at point source
70     u_pc=L.getValue(u)     u_pc=L.getValue(u)
71     print "u at point charge=",u_pc     print "u at point charge=",u_pc

u_pc_x = u_pc[0]
u_pc_y = u_pc[1]
u_pc_z = u_pc[2]

72     # open file to save displacement at point source     # open file to save displacement at point source
73     u_pc_data=FileWriter('./data/U_pc.out')     u_pc_data=FileWriter('./data/U_pc.out')
74     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))     u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
75
76     while t<tend:     while t<tend:
77         t+=h
78       # ... get current stress ....       # ... get current stress ....
80       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
81       # ... get new acceleration ....       # ... get new acceleration ....
82       mypde.setValue(X=-stress)                 amplitude=U0*2*exp(1)/alpha**2*(1-5*(t/alpha)**2+2*(t/alpha)**4)*exp(-(t/alpha)**2)
83         mypde.setValue(X=-stress, r=dunit*amplitude)
84       a=mypde.getSolution()       a=mypde.getSolution()
85       # ... get new displacement ...       # ... get new displacement ...
86       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
87       # ... shift displacements ....       # ... shift displacements ....
88       u_last=u       u_last=u
89       u=u_new       u=u_new
t+=h
90       n+=1       n+=1
91       print n,"-th time step t ",t       print n,"-th time step t ",t
92       u_pc=L.getValue(u)       u_pc=L.getValue(u)
93       print "u at point charge=",u_pc       print "u at point charge=",u_pc

u_pc_x=u_pc[0]
u_pc_y=u_pc[1]
u_pc_z=u_pc[2]

94       # save displacements at point source to file for t > 0       # save displacements at point source to file for t > 0
95       u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))       u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
96
97       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
98       if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,       if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
# Line 113  def wavePropagation(domain,h,tend,lam,mu Line 100  def wavePropagation(domain,h,tend,lam,mu
100
101     u_pc_data.close()     u_pc_data.close()
102
103  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
104    h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
105    print "time step size = ",h
106  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
107

Legend:
 Removed from v.2484 changed lines Added in v.2502