/[escript]/trunk/doc/examples/usersguide/wave.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2484 by gross, Mon Jun 22 04:22:19 2009 UTC revision 2502 by gross, Tue Jun 30 05:49:22 2009 UTC
# Line 19  __license__="""Licensed under the Open S Line 19  __license__="""Licensed under the Open S
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
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
    src_radius = 0.1*Lsup(domain.getSize())  
54     print "src_radius = ",src_radius     print "src_radius = ",src_radius
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    
58     mypde.setValue(D=kronecker*rho)     mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
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
63     u=U0*whereNegative(length(x-xc)-src_radius)*dunit     u=Vector(0.,Solution(domain))
64     u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit     u_last=Vector(0.,Solution(domain))
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 ....
79       g=grad(u)       g=grad(u)
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

  ViewVC Help
Powered by ViewVC 1.1.26