/[escript]/trunk/doc/cookbook/cblib/wavesolver2d.py
ViewVC logotype

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

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

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
57     src_radius = 0.1*Lsup(domain.getSize())     src_radius = 50#2*Lsup(domain.getSize())
58     print "src_radius = ",src_radius     print "src_radius = ",src_radius
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
67     u=U0*whereNegative(length(x-xc)-src_radius)*dunit     u=U0*whereNegative(length(x-xc)-src_radius)*dunit
68       print u
69     u_m1=U0*whereNegative(length(x-xc)-src_radius)*dunit     u_m1=U0*whereNegative(length(x-xc)-src_radius)*dunit
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
86         #g=grad(u)
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
93       g=grad(u)       g=grad(u)
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

  ViewVC Help
Powered by ViewVC 1.1.26