# Diff of /branches/stage3.0/doc/examples/usersguide/wave.py

revision 2583 by jfenwick, Fri Jul 31 05:39:36 2009 UTC revision 2584 by jfenwick, Wed Aug 5 02:44:51 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
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
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.2583 changed lines Added in v.2584