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

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

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

revision 3056 by ahallam, Sun Jul 4 21:52:36 2010 UTC revision 3057 by ahallam, Tue Jul 6 02:10:59 2010 UTC
# Line 29  __url__="https://launchpad.net/escript-f Line 29  __url__="https://launchpad.net/escript-f
29  #######################################################EXTERNAL MODULES  #######################################################EXTERNAL MODULES
30  from esys.escript import *  from esys.escript import *
31  from esys.finley import Rectangle  from esys.finley import Rectangle
 import sys  
32  import os  import os
33  # smoothing operator  # smoothing operator
34  from esys.escript.pdetools import Projector, Locator  from esys.escript.pdetools import Projector, Locator
# Line 39  import pylab as pl Line 38  import pylab as pl
38  import matplotlib.cm as cm  import matplotlib.cm as cm
39  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
40    
41    ########################################################MPI WORLD CHECK
42    if getMPISizeWorld() > 1:
43        import sys
44        print "This example will not run in an MPI world."
45        sys.exit(0)
46    
47  #################################################ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
48  # where to save output data  # where to save output data
49  savepath = "data/example08b"  savepath = "data/example08b"
# Line 46  mkDir(savepath) Line 51  mkDir(savepath)
51  #Geometric and material property related variables.  #Geometric and material property related variables.
52  mx = 1000. # model lenght  mx = 1000. # model lenght
53  my = 1000. # model width  my = 1000. # model width
54  ndx = 1000 # steps in x direction  ndx = 500 # steps in x direction
55  ndy = 1000 # steps in y direction  ndy = 500 # steps in y direction
56  xstep=mx/ndx # calculate the size of delta x  xstep=mx/ndx # calculate the size of delta x
57  ystep=abs(my/ndy) # calculate the size of delta y  ystep=abs(my/ndy) # calculate the size of delta y
58  lam=3.462e9 #lames constant  lam=3.462e9 #lames constant
# Line 55  mu=3.462e9  #bulk modulus Line 60  mu=3.462e9  #bulk modulus
60  rho=1154.   #density  rho=1154.   #density
61  # Time related variables.  # Time related variables.
62  tend=0.5    # end time  tend=0.5    # end time
63  h=0.0001    # time step  h=0.0005    # time step
64  # data recording times  # data recording times
65  rtime=0.0 # first time to record  rtime=0.0 # first time to record
66  rtime_inc=tend/50.0 # time increment to record  rtime_inc=tend/50.0 # time increment to record
# Line 74  dfeq=50 #Dominant Frequency Line 79  dfeq=50 #Dominant Frequency
79  a = 2.0 * (np.pi * dfeq)**2.0  a = 2.0 * (np.pi * dfeq)**2.0
80  t0 = 5.0 / (2.0 * np.pi * dfeq)  t0 = 5.0 / (2.0 * np.pi * dfeq)
81  srclength = 5. * t0  srclength = 5. * t0
82  ls = srclength/h  ls = int(srclength/h)
83  print 'source length',ls  print 'source length',ls
84  source=np.zeros(ls,'float') # source array  source=np.zeros(ls,'float') # source array
85  ampmax=0  ampmax=0
# Line 122  bleft=xstep*bn; bright=mx-(xstep*bn); bb Line 127  bleft=xstep*bn; bright=mx-(xstep*bn); bb
127  # locate these points in the domain  # locate these points in the domain
128  left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot  left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot
129    
130  tgamma=0.99   # decay value for exponential function  tgamma=0.85   # decay value for exponential function
131  def calc_gamma(G,npts):  def calc_gamma(G,npts):
132      func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))      func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))
133      return func      return func
# Line 178  t=0 # time counter Line 183  t=0 # time counter
183  ##############################################################ITERATION  ##############################################################ITERATION
184  while t<tend:  while t<tend:
185      # get current stress      # get current stress
186      g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))      g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))*abc
187      mypde.setValue(X=-stress) # set PDE values      mypde.setValue(X=-stress) # set PDE values
188      accel = mypde.getSolution() #get PDE solution for accelleration      accel = mypde.getSolution() #get PDE solution for accelleration
189      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement

Legend:
Removed from v.3056  
changed lines
  Added in v.3057

  ViewVC Help
Powered by ViewVC 1.1.26