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

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