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

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
# 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