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

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

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

revision 3068 by ahallam, Sun Jul 4 21:52:36 2010 UTC revision 3069 by ahallam, Wed Jul 21 03:24:48 2010 UTC
# Line 36  import numpy as np Line 36  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    from esys.finley import ReadMesh
40    
41  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
42  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
# Line 45  if getMPISizeWorld() > 1: Line 46  if getMPISizeWorld() > 1:
46    
47  #################################################ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
48  # where to save output data  # where to save output data
49  savepath = "data/example08b"  savepath = "data/example09"
50  mkDir(savepath)  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 = 500 # steps in x direction  mz=200.0
55  ndy = 500 # steps in y direction  step=5.0 # the element size
56  xstep=mx/ndx # calculate the size of delta x  ndx = int(mx/step) # steps in x direction
57  ystep=abs(my/ndy) # calculate the size of delta y  ndy = int(my/step) # steps in y direction
58  lam=3.462e9 #lames constant  ndz = int(mz/step)
59  mu=3.462e9  #bulk modulus  
60  rho=1154.   #density  vel2=1800.;   vel1=3000.
61    rho2=2300.;   rho1=3100. #density
62    mu2=(vel2**2.)*rho2/8.;  mu1=(vel1**2.)*rho1/8.  #bulk modulus
63    lam2=mu2*6.; lam1=mu1*6. #lames constant
64    
65  # Time related variables.  # Time related variables.
66  tend=0.5    # end time  tend=0.5    # end time
67  h=0.0005    # time step  h=0.0005    # time step
# Line 103  pl.plot(source) Line 108  pl.plot(source)
108  pl.savefig(os.path.join(savepath,'source.png'))  pl.savefig(os.path.join(savepath,'source.png'))
109    
110  # will introduce a spherical source at middle left of bottom face  # will introduce a spherical source at middle left of bottom face
111  xc=[mx/2,0]  xc=[mx/2,my/2,0.]
112    
113  ####################################################DOMAIN CONSTRUCTION  ####################################################DOMAIN CONSTRUCTION
114  domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) # create the domain  domain=ReadMesh(os.path.join(savepath,'example09m.fly')) # create the domain
115  x=domain.getX() # get the locations of the nodes in the domani  x=domain.getX() # get the locations of the nodes in the domain
116    
117    lam=Scalar(0,Function(domain))
118    lam.setTaggedValue("vintfa",lam1)
119    lam.setTaggedValue("vintfb",lam2)
120    mu=Scalar(0,Function(domain))
121    mu.setTaggedValue("vintfa",mu1)
122    mu.setTaggedValue("vintfb",mu2)
123    rho=Scalar(0,Function(domain))
124    rho.setTaggedValue("vintfa",rho1)
125    rho.setTaggedValue("vintfb",rho2)
126    
127  ##########################################################ESTABLISH PDE  ##########################################################ESTABLISH PDE
128  mypde=LinearPDE(domain) # create pde  mypde=LinearPDE(domain) # create pde
# Line 115  mypde.setSymmetryOn() # turn symmetry on Line 130  mypde.setSymmetryOn() # turn symmetry on
130  # turn lumping on for more efficient solving  # turn lumping on for more efficient solving
131  mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)  mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
132  kmat = kronecker(domain) # create the kronecker delta function of the domain  kmat = kronecker(domain) # create the kronecker delta function of the domain
133  mypde.setValue(D=kmat*rho) #set the general form value D  mypde.setValue(D=rho*kmat) #set the general form value D
134    
135    
136    
137  ##########################################################ESTABLISH ABC  ##########################################################ESTABLISH ABC
138  # Define where the boundary decay will be applied.  # Define where the boundary decay will be applied.
139  bn=50.  #bn=50.
140  bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn)  #bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn)
141  # btop=ystep*bn # don't apply to force boundary!!!  ## btop=ystep*bn # don't apply to force boundary!!!
142    
143  # locate these points in the domain  ## locate these points in the domain
144  left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot  #left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot
145    
146  tgamma=0.85   # decay value for exponential function  #tgamma=0.85   # decay value for exponential function
147  def calc_gamma(G,npts):  #def calc_gamma(G,npts):
148      func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))  #    func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))
149      return func  #    return func
150    
151  gleft  = calc_gamma(tgamma,bleft)  #gleft  = calc_gamma(tgamma,bleft)
152  gright = calc_gamma(tgamma,bleft)  #gright = calc_gamma(tgamma,bleft)
153  gbottom= calc_gamma(tgamma,ystep*bn)  #gbottom= calc_gamma(tgamma,ystep*bn)
154    
155  print 'gamma', gleft,gright,gbottom  #print 'gamma', gleft,gright,gbottom
156    
157  # calculate decay functions  ## calculate decay functions
158  def abc_bfunc(gamma,loc,x,G):  #def abc_bfunc(gamma,loc,x,G):
159      func=exp(-1.*(gamma*abs(loc-x))**2.)  #    func=exp(-1.*(gamma*abs(loc-x))**2.)
160      return func  #    return func
161    
162  fleft=abc_bfunc(gleft,bleft,x[0],tgamma)  #fleft=abc_bfunc(gleft,bleft,x[0],tgamma)
163  fright=abc_bfunc(gright,bright,x[0],tgamma)  #fright=abc_bfunc(gright,bright,x[0],tgamma)
164  fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma)  #fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma)
165  # apply these functions only where relevant  ## apply these functions only where relevant
166  abcleft=fleft*whereNegative(left)  #abcleft=fleft*whereNegative(left)
167  abcright=fright*wherePositive(right)  #abcright=fright*wherePositive(right)
168  abcbottom=fbottom*wherePositive(bottom)  #abcbottom=fbottom*wherePositive(bottom)
169  # make sure the inside of the abc is value 1  ## make sure the inside of the abc is value 1
170  abcleft=abcleft+whereZero(abcleft)  #abcleft=abcleft+whereZero(abcleft)
171  abcright=abcright+whereZero(abcright)  #abcright=abcright+whereZero(abcright)
172  abcbottom=abcbottom+whereZero(abcbottom)  #abcbottom=abcbottom+whereZero(abcbottom)
173  # multiply the conditions together to get a smooth result  ## multiply the conditions together to get a smooth result
174  abc=abcleft*abcright*abcbottom  #abc=abcleft*abcright*abcbottom
175    
176  #visualise the boundary function  #visualise the boundary function
177  abcT=abc.toListOfTuples()  #abcT=abc.toListOfTuples()
178  abcT=np.reshape(abcT,(ndx+1,ndy+1))  #abcT=np.reshape(abcT,(ndx+1,ndy+1))
179  pl.clf(); pl.imshow(abcT); pl.colorbar();  #pl.clf(); pl.imshow(abcT); pl.colorbar();
180  pl.savefig(os.path.join(savepath,"abc.png"))  #pl.savefig(os.path.join(savepath,"abc.png"))
181    
182    
183  ############################################FIRST TIME STEPS AND SOURCE  ############################################FIRST TIME STEPS AND SOURCE
184  # define small radius around point xc  # define small radius around point xc
185  src_length = 40; print "src_length = ",src_length  src_length = 40; print "src_length = ",src_length
186  # set initial values for first two time steps with source terms  # set initial values for first two time steps with source terms
187  y=source[0]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)  xb=FunctionOnBoundary(domain).getX()
188  src_dir=numpy.array([0.,-1.]) # defines direction of point source as down  y=source[0]*(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))
189    src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
190  y=y*src_dir  y=y*src_dir
191  mypde.setValue(y=y) #set the source as a function on the boundary  mypde.setValue(y=y) #set the source as a function on the boundary
192  # initial value of displacement at point source is constant (U0=0.01)  # initial value of displacement at point source is constant (U0=0.01)
193  # for first two time steps  # for first two time steps
194  u=[0.0,0.0]*whereNegative(x)  u=[0.0,0.0,0.0]*whereNegative(x)
195  u_m1=u  u_m1=u
196    
197  ####################################################ITERATION VARIABLES  ####################################################ITERATION VARIABLES
# Line 182  t=0 # time counter Line 200  t=0 # time counter
200  ##############################################################ITERATION  ##############################################################ITERATION
201  while t<tend:  while t<tend:
202      # get current stress      # get current stress
203      g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))*abc      g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
204      mypde.setValue(X=-stress) # set PDE values      mypde.setValue(X=-stress) # set PDE values
205      accel = mypde.getSolution() #get PDE solution for accelleration      accel = mypde.getSolution() #get PDE solution for accelleration
206      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
207      u_p1=u_p1*abc       # apply boundary conditions      u_p1=u_p1#*abc          # apply boundary conditions
208      u_m1=u; u=u_p1 # shift values by 1      u_m1=u; u=u_p1 # shift values by 1
209      # save current displacement, acceleration and pressure      # save current displacement, acceleration and pressure
210      if (t >= rtime):      if (t >= rtime):
211          saveVTK(os.path.join(savepath,"ex08b.%05d.vtu"%n),displacement=length(u),\          saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\
212                                      acceleration=length(accel),tensor=stress)                                      acceleration=length(accel),tensor=stress)
213          rtime=rtime+rtime_inc #increment data save time          rtime=rtime+rtime_inc #increment data save time
214      # increment loop values      # increment loop values

Legend:
Removed from v.3068  
changed lines
  Added in v.3069

  ViewVC Help
Powered by ViewVC 1.1.26