/[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 2549 by jfenwick, Mon Jul 20 06:43:47 2009 UTC revision 2580 by gross, Tue Aug 4 08:24:12 2009 UTC
# Line 24  from esys.escript.pdetools import Locato Line 24  from esys.escript.pdetools import Locato
24  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
25  from esys.finley import Brick  from esys.finley import Brick
26  from numpy import identity,zeros,ones  from numpy import identity,zeros,ones
27    import matplotlib.pyplot as plt
28    
29    
30  ne=32          # number of cells in x_0 and x_1 directions  ne=32          # number of cells in x_0 and x_1 directions
31  width=10000.  # length in x_0 and x_1 directions  width=10000.  # length in x_0 and x_1 directions
32  lam=3.462e9  lam=3.462e9
33  mu=3.462e9  mu=3.462e9
34  rho=1154.  rho=1154.
35  tend=20. # to ran a full simulation change tend to 60.  tend=10. # to ran a full simulation change tend to 60.
36  alpha=0.7  alpha=0.7
37  t0=3.  t0=3.
38    
39    
40  U0=1. # maximum displacement  U0=1. # maximum displacement
41  mkDir("data") # create directory data if it does not exist already.  mkDir("data") # create directory data if it does not exist already.
42    
43  def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):  def wavePropagation(domain,h,tend,lam,mu,rho, xc, src_radius, U0):
44       # lists to collect displacement at point source
45       ts, u_pc0,u_pc1,u_pc2=[], [], [], []
46     x=domain.getX()     x=domain.getX()
47     # ... open new PDE ...     # ... open new PDE ...
48     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
# Line 59  def wavePropagation(domain,h,tend,lam,mu Line 64  def wavePropagation(domain,h,tend,lam,mu
64     # find potential at point source     # find potential at point source
65     u_pc=L.getValue(u)     u_pc=L.getValue(u)
66     print "u at point charge=",u_pc     print "u at point charge=",u_pc
67     # open file to save displacement at point source     ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
    u_pc_data=FileWriter('./data/U_pc.out')  
    u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))  
68    
69     while t<tend:     while t<tend:
70       t+=h       t+=h
# Line 81  def wavePropagation(domain,h,tend,lam,mu Line 84  def wavePropagation(domain,h,tend,lam,mu
84       print n,"-th time step t ",t       print n,"-th time step t ",t
85       u_pc=L.getValue(u)       u_pc=L.getValue(u)
86       print "u at point charge=",u_pc       print "u at point charge=",u_pc
87       # save displacements at point source to file for t > 0       ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
      u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))  
88    
89       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
90       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,
91       displacement = length(u), tensor = stress, Ux = u[0] )       displacement = length(u), tensor = stress, Ux = u[0] )
92       return ts, u_pc0,u_pc1,u_pc2
93    
94     u_pc_data.close()  #
95      # create domain:
96    #
97  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/ne)
98    #
99    #  sety time step size:
100    #
101  h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())  h=inf(1./5.)*inf(sqrt(rho/(lam+2*mu))*mydomain.getSize())
102  print "time step size = ",h  print "time step size = ",h
103    #
104  #  spherical source at middle of bottom face  #  spherical source at middle of bottom face
105    #
106  xc=[width/2.,width/2.,0.]  xc=[width/2.,width/2.,0.]
107  # define small radius around point xc  # define small radius around point xc
108  src_radius = 0.03*width  src_radius = 0.03*width
109  print "src_radius = ",src_radius  print "src_radius = ",src_radius
110  wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)  #
111    # run it
112    #
113    ts, u_pc0,u_pc1,u_pc2 = wavePropagation(mydomain,h,tend,lam,mu,rho, xc, src_radius, U0)
114    #
115    # create a plot:
116    #
117    if getMPIRankWorld() == 0:
118        plt.title("Displacement at Point Source")
119        plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
120        plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
121        plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
122        plt.xlabel('time')
123        plt.ylabel('displacement')
124        plt.legend()
125        plt.savefig('u_pc.png', format='png')
126    # or save displacement
127    u_pc_data=FileWriter('./data/U_pc.out')
128    for i in xrange(len(ts)) :
129        u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
130    u_pc_data.close()
131    

Legend:
Removed from v.2549  
changed lines
  Added in v.2580

  ViewVC Help
Powered by ViewVC 1.1.26