/[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 3089 by ahallam, Mon Aug 9 07:20:58 2010 UTC
# Line 49  if getMPISizeWorld() > 1: Line 49  if getMPISizeWorld() > 1:
49  savepath = "data/example09"  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 = 200. # model lenght
53  my = 1000. # model width  my = 200. # model width
54  mz=200.0  mz=100.0
55  step=5.0 # the element size  step=10.0 # the element size
56  ndx = int(mx/step) # steps in x direction  ndx = int(mx/step) # steps in x direction
57  ndy = int(my/step) # steps in y direction  ndy = int(my/step) # steps in y direction
58  ndz = int(mz/step)  ndz = int(mz/step)
59    
60  vel2=1800.;   vel1=3000.  vel2=1800.;   vel1=3000.
61  rho2=2300.;   rho1=3100. #density  rho2=2300.;   rho1=3100. #density
62  mu2=(vel2**2.)*rho2/8.;  mu1=(vel1**2.)*rho1/8.  #bulk modulus  mu2=vel2**2.*rho2/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus
63  lam2=mu2*6.; lam1=mu1*6. #lames constant  lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant
64    
65  # Time related variables.  # Time related variables.
66  tend=0.5    # end time  tend=0.1    # end time
67  h=0.0005    # time step  h=0.0001    # time step
68  # data recording times  # data recording times
69  rtime=0.0 # first time to record  rtime=0.0 # first time to record
70  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 92  for it in range(0,ls):
92      tt = t-t0      tt = t-t0
93      dum1 = np.exp(-a * tt * tt)      dum1 = np.exp(-a * tt * tt)
94      source[it] = -2. * a * tt * dum1      source[it] = -2. * a * tt * dum1
 #   source[it] = exp(-a * tt * tt)    !gaussian  
95      if (abs(source[it]) > ampmax):      if (abs(source[it]) > ampmax):
96          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)  
97      time[t]=t*h      time[t]=t*h
98  #tdecay=decay1*decay2*U0  #tdecay=decay1*decay2*U0
99  #decay1=decay1*U0; decay2=decay2*U0  #decay1=decay1*U0; decay2=decay2*U0
100  pl.clf();  #pl.clf();
101  pl.plot(source)  #pl.plot(source)
102  #pl.plot(time,decay1);pl.plot(time,decay2);  #pl.plot(time,decay1);pl.plot(time,decay2);
103  #pl.plot(time,tdecay)  #pl.plot(time,tdecay)
104  pl.savefig(os.path.join(savepath,'source.png'))  #pl.savefig(os.path.join(savepath,'source.png'))
105    
106  # will introduce a spherical source at middle left of bottom face  # will introduce a spherical source at middle left of bottom face
107  xc=[mx/2,my/2,0.]  xc=[50,50,0]
108    
109  ####################################################DOMAIN CONSTRUCTION  ####################################################DOMAIN CONSTRUCTION
110  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 128  mypde.setSymmetryOn() # turn symmetry on
128  kmat = kronecker(domain) # create the kronecker delta function of the domain  kmat = kronecker(domain) # create the kronecker delta function of the domain
129  mypde.setValue(D=rho*kmat) #set the general form value D  mypde.setValue(D=rho*kmat) #set the general form value D
130    
   
   
 ##########################################################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"))  
   
   
131  ############################################FIRST TIME STEPS AND SOURCE  ############################################FIRST TIME STEPS AND SOURCE
132  # define small radius around point xc  # define small radius around point xc
133  src_length = 40; print "src_length = ",src_length  src_length = 20; print "src_length = ",src_length
134  # set initial values for first two time steps with source terms  # set initial values for first two time steps with source terms
135  xb=FunctionOnBoundary(domain).getX()  xb=FunctionOnBoundary(domain).getX()
136  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))
137    y=Vector(0.0,FunctionOnBoundary(domain))
138    
139  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
140  y=y*src_dir  
141    #sy=sy*src_dir
142    #sy.setTaggedValue("stop")
143    y.setTaggedValue("stop",src_dir*source[0])#)*(cos(length(xb-xc)*3.1415/src_length)+1))
144  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
145  # initial value of displacement at point source is constant (U0=0.01)  # initial value of displacement at point source is constant (U0=0.01)
146  # for first two time steps  # for first two time steps
147  u=[0.0,0.0,0.0]*whereNegative(x)  u=[0.0,0.0,0.0]*wherePositive(x)
148  u_m1=u  u_m1=u
149    
150  ####################################################ITERATION VARIABLES  ####################################################ITERATION VARIABLES
# Line 214  while t<tend: Line 167  while t<tend:
167      # increment loop values      # increment loop values
168      t=t+h; n=n+1      t=t+h; n=n+1
169      if (n < ls):      if (n < ls):
170          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)
171          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
172      print n,"-th time step t ",t      print n,"-th time step t ",t

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

  ViewVC Help
Powered by ViewVC 1.1.26