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

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

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

trunk/doc/examples/cookbook/example09.py revision 3075 by ahallam, Wed Jul 28 02:51:20 2010 UTC trunk/doc/examples/cookbook/example09a.py revision 3346 by caltinay, Fri Nov 12 01:19:02 2010 UTC
# Line 33  import os Line 33  import os
33  from esys.escript.pdetools import Projector, Locator  from esys.escript.pdetools import Projector, Locator
34  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
35  import numpy as np  import numpy as np
36    import matplotlib
37    matplotlib.use('agg') #It's just here for automated testing
38    
39  import pylab as pl  import pylab as pl
40  import matplotlib.cm as cm  import matplotlib.cm as cm
41  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
42  from esys.finley import ReadMesh  from esys.finley import ReadMesh
43    from esys.weipa import saveVTK
44    
45  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
46  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
# Line 49  if getMPISizeWorld() > 1: Line 53  if getMPISizeWorld() > 1:
53  savepath = "data/example09"  savepath = "data/example09"
54  mkDir(savepath)  mkDir(savepath)
55  #Geometric and material property related variables.  #Geometric and material property related variables.
56  mx = 1000. # model lenght  mx = 200. # model lenght
57  my = 1000. # model width  my = 200. # model width
58  mz=200.0  mz=100.0
59  step=5.0 # the element size  step=10.0 # the element size
60  ndx = int(mx/step) # steps in x direction  ndx = int(mx/step) # steps in x direction
61  ndy = int(my/step) # steps in y direction  ndy = int(my/step) # steps in y direction
62  ndz = int(mz/step)  ndz = int(mz/step)
63    
64  vel2=1800.;   vel1=3000.  vel2=1800.;   vel1=3000.
65  rho2=2300.;   rho1=3100. #density  rho2=2300.;   rho1=3100. #density
66  mu2=(vel2**2.)*rho2/8.;  mu1=(vel1**2.)*rho1/8.  #bulk modulus  mu2=vel2**2.*rho2/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus
67  lam2=mu2*6.; lam1=mu1*6. #lames constant  lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant
68    
69  # Time related variables.  # Time related variables.
70  tend=0.5    # end time  testing=True
71  h=0.0005    # time step  if testing:
72        print 'The testing end time is curerntly sellected this severely limits the number of time iterations.'
73        print "Try changing testing to False for more iterations."
74        tend=0.001
75    else:
76        tend=0.1    # end time
77    
78    h=0.0001    # time step
79  # data recording times  # data recording times
80  rtime=0.0 # first time to record  rtime=0.0 # first time to record
81  rtime_inc=tend/50.0 # time increment to record  rtime_inc=tend/50.0 # time increment to record
# Line 92  for it in range(0,ls): Line 103  for it in range(0,ls):
103      tt = t-t0      tt = t-t0
104      dum1 = np.exp(-a * tt * tt)      dum1 = np.exp(-a * tt * tt)
105      source[it] = -2. * a * tt * dum1      source[it] = -2. * a * tt * dum1
 #   source[it] = exp(-a * tt * tt)    !gaussian  
106      if (abs(source[it]) > ampmax):      if (abs(source[it]) > ampmax):
107          ampmax = abs(source[it])          ampmax = abs(source[it])
     #source[t]=np.exp(g*t)*U0*np.sin(2.*np.pi*t/(0.75*ls))*(np.exp(-.1*g*t)-1)  
     #decay1[t]=np.exp(g*t)  
     #decay2[t]=(np.exp(-.1*g*t)-1)  
108      time[t]=t*h      time[t]=t*h
109  #tdecay=decay1*decay2*U0  #tdecay=decay1*decay2*U0
110  #decay1=decay1*U0; decay2=decay2*U0  #decay1=decay1*U0; decay2=decay2*U0
111  pl.clf();  #pl.clf();
112  pl.plot(source)  #pl.plot(source)
113  #pl.plot(time,decay1);pl.plot(time,decay2);  #pl.plot(time,decay1);pl.plot(time,decay2);
114  #pl.plot(time,tdecay)  #pl.plot(time,tdecay)
115  pl.savefig(os.path.join(savepath,'source.png'))  #pl.savefig(os.path.join(savepath,'source.png'))
116    
117  # will introduce a spherical source at middle left of bottom face  # will introduce a spherical source at middle left of bottom face
118  xc=[mx/2,my/2,0.]  xc=[50,50,0]
119    
120  ####################################################DOMAIN CONSTRUCTION  ####################################################DOMAIN CONSTRUCTION
121  domain=ReadMesh(os.path.join(savepath,'example09m.fly')) # create the domain  domain=ReadMesh(os.path.join(savepath,'example09m.fly')) # create the domain
# Line 132  mypde.setSymmetryOn() # turn symmetry on Line 139  mypde.setSymmetryOn() # turn symmetry on
139  kmat = kronecker(domain) # create the kronecker delta function of the domain  kmat = kronecker(domain) # create the kronecker delta function of the domain
140  mypde.setValue(D=rho*kmat) #set the general form value D  mypde.setValue(D=rho*kmat) #set the general form value D
141    
   
   
 ##########################################################ESTABLISH ABC  
 # Define where the boundary decay will be applied.  
 #bn=50.  
 #bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn)  
 ## btop=ystep*bn # don't apply to force boundary!!!  
   
 ## locate these points in the domain  
 #left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot  
   
 #tgamma=0.85   # decay value for exponential function  
 #def calc_gamma(G,npts):  
 #    func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))  
 #    return func  
   
 #gleft  = calc_gamma(tgamma,bleft)  
 #gright = calc_gamma(tgamma,bleft)  
 #gbottom= calc_gamma(tgamma,ystep*bn)  
   
 #print 'gamma', gleft,gright,gbottom  
   
 ## calculate decay functions  
 #def abc_bfunc(gamma,loc,x,G):  
 #    func=exp(-1.*(gamma*abs(loc-x))**2.)  
 #    return func  
   
 #fleft=abc_bfunc(gleft,bleft,x[0],tgamma)  
 #fright=abc_bfunc(gright,bright,x[0],tgamma)  
 #fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma)  
 ## apply these functions only where relevant  
 #abcleft=fleft*whereNegative(left)  
 #abcright=fright*wherePositive(right)  
 #abcbottom=fbottom*wherePositive(bottom)  
 ## make sure the inside of the abc is value 1  
 #abcleft=abcleft+whereZero(abcleft)  
 #abcright=abcright+whereZero(abcright)  
 #abcbottom=abcbottom+whereZero(abcbottom)  
 ## multiply the conditions together to get a smooth result  
 #abc=abcleft*abcright*abcbottom  
   
 #visualise the boundary function  
 #abcT=abc.toListOfTuples()  
 #abcT=np.reshape(abcT,(ndx+1,ndy+1))  
 #pl.clf(); pl.imshow(abcT); pl.colorbar();  
 #pl.savefig(os.path.join(savepath,"abc.png"))  
   
   
142  ############################################FIRST TIME STEPS AND SOURCE  ############################################FIRST TIME STEPS AND SOURCE
143  # define small radius around point xc  # define small radius around point xc
144  src_length = 40; print "src_length = ",src_length  src_length = 20; print "src_length = ",src_length
145  # set initial values for first two time steps with source terms  # set initial values for first two time steps with source terms
146  xb=FunctionOnBoundary(domain).getX()  xb=FunctionOnBoundary(domain).getX()
147  y=source[0]*(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))  #sy=source[0]*(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))
148    y=Vector(0.0,FunctionOnBoundary(domain))
149    
150  src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down  src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
151  y=y*src_dir  
152    #sy=sy*src_dir
153    #sy.setTaggedValue("stop")
154    y.setTaggedValue("stop",src_dir*source[0])#)*(cos(length(xb-xc)*3.1415/src_length)+1))
155  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
156  # initial value of displacement at point source is constant (U0=0.01)  # initial value of displacement at point source is constant (U0=0.01)
157  # for first two time steps  # for first two time steps
158  u=[0.0,0.0,0.0]*whereNegative(x)  u=[0.0,0.0,0.0]*wherePositive(x)
159  u_m1=u  u_m1=u
160    
161  ####################################################ITERATION VARIABLES  ####################################################ITERATION VARIABLES
# Line 214  while t<tend: Line 178  while t<tend:
178      # increment loop values      # increment loop values
179      t=t+h; n=n+1      t=t+h; n=n+1
180      if (n < ls):      if (n < ls):
181          y=source[n]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)          y.setTaggedValue("stop",source[n]*src_dir)
182          y=y*src_dir; 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
183      print n,"-th time step t ",t      print n,"-th time step t ",t

Legend:
Removed from v.3075  
changed lines
  Added in v.3346

  ViewVC Help
Powered by ViewVC 1.1.26