/[escript]/trunk/doc/examples/cookbook/example07a.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/example07a.py

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

revision 3004 by ahallam, Wed Apr 14 05:27:30 2010 UTC revision 3025 by ahallam, Thu May 6 01:20:46 2010 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    
22    ############################################################FILE HEADER
23    # example07a.py
24  # Antony Hallam  # Antony Hallam
25  # Acoustic Wave Equation Simulation  # Acoustic Wave Equation Simulation using displacement solution
26    
27  # Importing all the necessary modules required.  #######################################################EXTERNAL MODULES
28  from esys.escript import *  from esys.escript import *
29  from esys.finley import Rectangle  from esys.finley import Rectangle
30  import sys  import sys
31  import os  import os
32  # smoothing operator  # smoothing operator
33  from esys.escript.pdetools import Projector  from esys.escript.pdetools import Projector, Locator
34    from esys.escript.unitsSI import *
35  import numpy as np  import numpy as np
36  import pylab as pl  import pylab as pl
37  import matplotlib.cm as cm  import matplotlib.cm as cm
38  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
39    
40  # Establish a save path.  ########################################################MPI WORLD CHECK
41    if getMPISizeWorld() > 1:
42        import sys
43        print "This example will not run in an MPI world."
44        sys.exit(0)
45    
46    #################################################ESTABLISHING VARIABLES
47    # where to save output data
48  savepath = "data/example07"  savepath = "data/example07"
49  mkDir(savepath)  mkDir(savepath)
   
50  #Geometric and material property related variables.  #Geometric and material property related variables.
51  mx = 1000. # model lenght  mx = 1000. # model lenght
52  my = 1000. # model width  my = 1000. # model width
53  ndx = 400 # steps in x direction  ndx = 400 # steps in x direction
54  ndy = 400 # steps in y direction  ndy = 400 # steps in y direction
55    xstep=mx/ndx # calculate the size of delta x
56    ystep=my/ndy # calculate the size of delta y
57    
58  xstep=mx/ndx  c=380.0*m/sec  # velocity of sound in air
59  ystep=my/ndy  csq=c*c #square of c
   
 c=380.0  
 csq=c*c  
60  # Time related variables.  # Time related variables.
61  tend=1.5    #end time  tend=1.5    # end time
62  # timestep  h=0.001     # time step
63  h=0.001  # data recording times
64    rtime=0.0 # first time to record
65    rtime_inc=tend/20.0 # time increment to record
66  #Check to make sure number of time steps is not too large.  #Check to make sure number of time steps is not too large.
67  print "Time step size= ",h, "Expected number of outputs= ",tend/h  print "Time step size= ",h, "Expected number of outputs= ",tend/h
68    
69  U0=0.01 # amplitude of point source  U0=0.005 # amplitude of point source
70  xc=[500,500] #location of point source  # want a spherical source in the middle of area
71    xc=[500,500] # with reference to mx,my this is the source location
72    
73    ####################################################DOMAIN CONSTRUCTION
74    mydomain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) # create the domain
75    x=mydomain.getX() # get the node locations of the domain
76    
77    ##########################################################ESTABLISH PDE
78    mypde=LinearPDE(mydomain) # create pde
79    mypde.setSymmetryOn() # turn symmetry on
80    mypde.setValue(D=1.) # set the value of D in the general form to 1.
81    
82  mydomain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  ############################################FIRST TIME STEPS AND SOURCE
 x=mydomain.getX()  
 # ... open new PDE ...  
 mypde=LinearPDE(mydomain)  
 mypde.setSymmetryOn()  
 mypde.setValue(D=1.)  
83  # define small radius around point xc  # define small radius around point xc
84  src_radius = 30  src_radius = 30
85  print "src_radius = ",src_radius  print "src_radius = ",src_radius
86    # set initial values for first two time steps with source terms
 # ... set initial values ....  
 n=0  
 # for first two time steps  
87  u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)  u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)
88  u_m1=u  u_m1=u
 t=0  
   
89  #plot source shape  #plot source shape
90  uT=np.array(u.toListOfTuples())  cut_loc=[] #where the cross section of the source along x will be
91  uT=np.reshape(uT,(ndx+1,ndy+1))  src_cut=[] #where the cross section of the source will be
92  source_line=uT[ndx/2,:]  # create locations for source cross section
93  pl.plot(source_line)  for i in range(ndx/2-ndx/10,ndx/2+ndx/10):
94  pl.plot(source_line,'ro')      cut_loc.append(xstep*i)
95  pl.axis([70,130,0,0.05])      src_cut.append([xstep*i,xc[1]])
96    # locate the nearest nodes to the points in src_cut
97    src=Locator(mydomain,src_cut)
98    src_cut=src.getValue(u) #retrieve the values from the nodes
99    # plot the x locations vs value and save the figure
100    pl.plot(cut_loc,src_cut)
101    pl.axis([xc[0]-src_radius*3,xc[0]+src_radius*3,0.,2.*U0])
102  pl.savefig(os.path.join(savepath,"source_line.png"))  pl.savefig(os.path.join(savepath,"source_line.png"))
103    
104    ####################################################ITERATION VARIABLES
105    n=0 # iteration counter
106    t=0 # time counter
107    ##############################################################ITERATION
108  while t<tend:  while t<tend:
109      # get current pressure      g=grad(u); pres=csq*h*h*g # get current pressure
110      g=grad(u)      mypde.setValue(X=-pres,Y=(2.*u-u_m1)) # set values in pde
111      pres=csq*h*h*g      u_p1 = mypde.getSolution() # get the new displacement
112      # set values and calculate solution      u_m1=u; u=u_p1 # shift values back one time step for next iteration
113      mypde.setValue(X=-pres,Y=(2.*u-u_m1))          # save current displacement, acceleration and pressure
114      u_p1 = mypde.getSolution()      if (t >= rtime):
115      # shift displacements          saveVTK(os.path.join(savepath,"ex07b.%i.vtu"%n),displacement=length(u),tensor=pres)
116      u_m1=u; u=u_p1          rtime=rtime+rtime_inc #increment data save time
117      # iteration increments      # increment loop values
118      t+=h; n+=1      t=t+h; n=n+1
119      print n,"-th time step t ",t      print n,"-th time step t ",t
     # ... save current acceleration in units of gravity and displacements  
     saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=pres)  

Legend:
Removed from v.3004  
changed lines
  Added in v.3025

  ViewVC Help
Powered by ViewVC 1.1.26