# Diff of /trunk/doc/cookbook/cblib/wavesolver2d.py

revision 2433 by ahallam, Fri May 15 04:29:01 2009 UTC revision 2434 by ahallam, Thu May 21 07:59:53 2009 UTC
# Line 29  from esys.finley import Rectangle Line 29  from esys.finley import Rectangle
29  from numarray import identity,zeros,ones  from numarray import identity,zeros,ones
30  import os  import os
31
32    ########################################################
33    # subroutine: wavesolver2d
34    # Can solve a generic 2D wave propagation problem with a
35    # point source in a homogeneous medium.
36    # Arguments:
37    #   domain  : domain to solve over
38    #   h       : time step
39    #   tend    : end time
40    #   lam, mu : lames constants for domain
41    #   rho : density of domain
42    #   U0  : magnitude of source
43    #   xc  : source location in domain (Vector)
44    #   savepath: where to output the data files
45    ########################################################
46  def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):  def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):
47     x=domain.getX()     x=domain.getX()
48     # ... open new PDE ...     # ... open new PDE ...
49     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
50     mypde.setSolverMethod(LinearPDE.LUMPING)    # mypde.setSolverMethod(LinearPDE.LUMPING)
51     kmat = kronecker(mypde.getDim())     mypde.setSymmetryOn()
52       kmat = kronecker(domain)
53       mypde.setValue(D=kmat*rho)
54
55     # define small radius around point xc     # define small radius around point xc
56     # Lsup(x) returns the maximum value of the argument x     # Lsup(x) returns the maximum value of the argument x
59
60     dunit=numarray.array([1.,0.]) # defines direction of point source     dunit=numarray.array([1.,0.]) # defines direction of point source
61
62     mypde.setValue(D=kmat*rho)
63     # ... set initial values ....     # ... set initial values ....
64     n=0     n=0
65     # initial value of displacement at point source is constant (U0=0.01)     # initial value of displacement at point source is constant (U0=0.01)
66     # for first two time steps     # for first two time steps
68       print u
70     t=0     t=0
71
# Line 65  def wavesolver2d(domain,h,tend,lam,mu,rh Line 82  def wavesolver2d(domain,h,tend,lam,mu,rh
82     u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')     u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
83     u_pc_data.write("%f %f %f\n"%(t,u_pc_x,u_pc_y))     u_pc_data.write("%f %f %f\n"%(t,u_pc_x,u_pc_y))
84
85         #NEW WAY
87         #stress=lam*trace(g)*kmat+mu*(g+transpose(g))
88
89     while t<tend:     while t<tend:
90       # ... get current stress ....       # ... get current stress ....
91
92    ##OLD WAY
94       stress=lam*trace(g)*kmat+mu*(g+transpose(g))       stress=lam*trace(g)*kmat+mu*(g+transpose(g))
95       # ... get new acceleration ....       ### ... get new acceleration ....
96       mypde.setValue(X=-stress)                 mypde.setValue(X=-stress)
97       a=mypde.getSolution()       a=mypde.getSolution()
98       # ... get new displacement ...       ### ... get new displacement ...
99       u_p1=2*u-u_last+h**2*a       u_p1=2*u-u_m1+h*h*a
100    ###NEW WAY
101         #mypde.setValue(X=-stress*h**2)
102         #mypde.setValue(Y=rho*(2*u-u_m1))
103         #u_p1 = mypde.getSolution()
104
105       # ... shift displacements ....       # ... shift displacements ....
106       u_m1=u       u_m1=u
107       u=u_p1       u=u_p1
108         #stress =
109       t+=h       t+=h
110       n+=1       n+=1
111       print n,"-th time step t ",t       print n,"-th time step t ",t
112       u_pc=L.getValue(u)       u_pc=L.getValue(u)
113       print "u at point charge=",u_pc  #     print "u at point charge=",u_pc
114
115       u_pc_x=u_pc[0]       u_pc_x=u_pc[0]
116       u_pc_y=u_pc[1]       u_pc_y=u_pc[1]
# Line 90  def wavesolver2d(domain,h,tend,lam,mu,rh Line 119  def wavesolver2d(domain,h,tend,lam,mu,rh
119       u_pc_data.write("%f %f %f\n"%(t,u_pc_x,u_pc_y))       u_pc_data.write("%f %f %f\n"%(t,u_pc_x,u_pc_y))
120
121       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
122       if n==1 or n%10==0: saveVTK(os.path.join(savepath,"usoln.%i.vtu"%(n/10)),acceleration=length(a)/9.81,       #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
123       displacement = length(u), tensor = stress, Ux = u[0] )       #displacement = length(u), tensor = stress, Ux = u[0] )
124         saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
125     u_pc_data.close()     u_pc_data.close()

Legend:
 Removed from v.2433 changed lines Added in v.2434