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

revision 5287 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 1  Line 1
1  from __future__ import division  from __future__ import division, print_function
from __future__ import print_function
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2009-2014 by University of Queensland  # Copyright (c) 2009-2014 by University of Queensland
32  import matplotlib  import matplotlib
33  matplotlib.use('agg') #It's just here for automated testing  matplotlib.use('agg') #It's just here for automated testing
34  from esys.escript import *  from esys.escript import *
from esys.finley import Rectangle
35  from esys.weipa import saveVTK  from esys.weipa import saveVTK
36  import os  import os
37  # smoothing operator  # smoothing operator
# Line 44  import numpy as np Line 42  import numpy as np
42  import pylab as pl  import pylab as pl
43  import matplotlib.cm as cm  import matplotlib.cm as cm
44  from esys.escript.linearPDEs import LinearPDE, SolverOptions  from esys.escript.linearPDEs import LinearPDE, SolverOptions
45    try:
46        # This imports the rectangle domain function
47        from esys.finley import Rectangle
48        HAVE_FINLEY = True
49    except ImportError:
50        print("Finley module not available")
51        HAVE_FINLEY = False
52  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
53  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
54      import sys      import sys
55      print("This example will not run in an MPI world.")      print("This example will not run in an MPI world.")
56      sys.exit(0)      sys.exit(0)
57
58  #################################################ESTABLISHING VARIABLES  if HAVE_FINLEY:
59  # where to save output data      #################################################ESTABLISHING VARIABLES
60  savepath = "data/example08b"      # where to save output data
61  mkDir(savepath)      savepath = "data/example08b"
62  #Geometric and material property related variables.      mkDir(savepath)
63  mx = 1000. # model lenght      #Geometric and material property related variables.
64  my = 1000. # model width      mx = 1000. # model lenght
65  ndx = 300 # steps in x direction      my = 1000. # model width
66  ndy = 300 # steps in y direction      ndx = 300 # steps in x direction
67  xstep=mx/ndx # calculate the size of delta x      ndy = 300 # steps in y direction
68  ystep=abs(my/ndy) # calculate the size of delta y      xstep=mx/ndx # calculate the size of delta x
69  lam=3.462e9 #lames constant      ystep=abs(my/ndy) # calculate the size of delta y
70  mu=3.462e9  #bulk modulus      lam=3.462e9 #lames constant
71  rho=1154.   #density      mu=3.462e9  #bulk modulus
72  # Time related variables.      rho=1154.   #density
73  testing=True      # Time related variables.
74  if testing:      testing=True
75      print('The testing end time is currently selected. This severely limits the number of time iterations.')      if testing:
76      print("Try changing testing to False for more iterations.")          print('The testing end time is currently selected. This severely limits the number of time iterations.')
77      tend=0.001          print("Try changing testing to False for more iterations.")
78  else:          tend=0.001
79      tend=0.5    # end time      else:
80            tend=0.5    # end time
81  h=0.0001    # time step
82  # data recording times      h=0.0001    # time step
83  rtime=0.0 # first time to record      # data recording times
84  rtime_inc=tend/50.0 # time increment to record      rtime=0.0 # first time to record
85  #Check to make sure number of time steps is not too large.      rtime_inc=tend/50.0 # time increment to record
86  print("Time step size= ",h, "Expected number of outputs= ",tend/h)      #Check to make sure number of time steps is not too large.
87        print("Time step size= ",h, "Expected number of outputs= ",tend/h)
88  U0=0.1 # amplitude of point source
89  ls=500   # length of the source      U0=0.1 # amplitude of point source
90  source=np.zeros(ls,'float') # source array      ls=500   # length of the source
91  decay1=np.zeros(ls,'float') # decay curve one      source=np.zeros(ls,'float') # source array
92  decay2=np.zeros(ls,'float') # decay curve two      decay1=np.zeros(ls,'float') # decay curve one
93  time=np.zeros(ls,'float')   # time values      decay2=np.zeros(ls,'float') # decay curve two
94  g=np.log(0.01)/ls      time=np.zeros(ls,'float')   # time values
95        g=np.log(0.01)/ls
96  dfeq=50 #Dominant Frequency
97  a = 2.0 * (np.pi * dfeq)**2.0      dfeq=50 #Dominant Frequency
98  t0 = 5.0 / (2.0 * np.pi * dfeq)      a = 2.0 * (np.pi * dfeq)**2.0
99  srclength = 5. * t0      t0 = 5.0 / (2.0 * np.pi * dfeq)
100  ls = int(srclength/h)      srclength = 5. * t0
101  print('source length',ls)      ls = int(srclength/h)
102  source=np.zeros(ls,'float') # source array      print('source length',ls)
103  ampmax=0      source=np.zeros(ls,'float') # source array
104  for it in range(0,ls):      ampmax=0
105      t = it*h      for it in range(0,ls):
106      tt = t-t0          t = it*h
107      dum1 = np.exp(-a * tt * tt)          tt = t-t0
108      source[it] = -2. * a * tt * dum1          dum1 = np.exp(-a * tt * tt)
109  #   source[it] = exp(-a * tt * tt)    !gaussian          source[it] = -2. * a * tt * dum1
110      if (abs(source[it]) > ampmax):      #   source[it] = exp(-a * tt * tt)    !gaussian
111          ampmax = abs(source[it])          if (abs(source[it]) > ampmax):
112      #source[t]=np.exp(g*t)*U0*np.sin(2.*np.pi*t/(0.75*ls))*(np.exp(-.1*g*t)-1)              ampmax = abs(source[it])
113      #decay1[t]=np.exp(g*t)          #source[t]=np.exp(g*t)*U0*np.sin(2.*np.pi*t/(0.75*ls))*(np.exp(-.1*g*t)-1)
114      #decay2[t]=(np.exp(-.1*g*t)-1)          #decay1[t]=np.exp(g*t)
115      time[t]=t*h          #decay2[t]=(np.exp(-.1*g*t)-1)
116  #tdecay=decay1*decay2*U0          time[t]=t*h
117  #decay1=decay1*U0; decay2=decay2*U0      #tdecay=decay1*decay2*U0
118  pl.clf();      #decay1=decay1*U0; decay2=decay2*U0
119  pl.plot(source)      pl.clf();
120  #pl.plot(time,decay1);pl.plot(time,decay2);      pl.plot(source)
121  #pl.plot(time,tdecay)      #pl.plot(time,decay1);pl.plot(time,decay2);
122  pl.savefig(os.path.join(savepath,'source.png'))      #pl.plot(time,tdecay)
123        pl.savefig(os.path.join(savepath,'source.png'))
124  # will introduce a spherical source at middle left of bottom face
125  xc=[mx/2,0]      # will introduce a spherical source at middle left of bottom face
126        xc=[mx/2,0]
127  ####################################################DOMAIN CONSTRUCTION
128  domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain      ####################################################DOMAIN CONSTRUCTION
129  x=domain.getX() # get the locations of the nodes in the domani      domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain
130        x=domain.getX() # get the locations of the nodes in the domani
131  ##########################################################ESTABLISH PDE
132  mypde=LinearPDE(domain) # create pde      ##########################################################ESTABLISH PDE
133  mypde.setSymmetryOn() # turn symmetry on      mypde=LinearPDE(domain) # create pde
134  # turn lumping on for more efficient solving      mypde.setSymmetryOn() # turn symmetry on
135  mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)      # turn lumping on for more efficient solving
136  kmat = kronecker(domain) # create the kronecker delta function of the domain      mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
137  mypde.setValue(D=kmat*rho) #set the general form value D      kmat = kronecker(domain) # create the kronecker delta function of the domain
138        mypde.setValue(D=kmat*rho) #set the general form value D
139  ##########################################################ESTABLISH ABC
140  # Define where the boundary decay will be applied.      ##########################################################ESTABLISH ABC
141  bn=50.      # Define where the boundary decay will be applied.
142  bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn)      bn=50.
143  # btop=ystep*bn # don't apply to force boundary!!!      bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn)
144        # btop=ystep*bn # don't apply to force boundary!!!
145  # locate these points in the domain
146  left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot      # locate these points in the domain
147        left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot
148  tgamma=0.85   # decay value for exponential function
149  def calc_gamma(G,npts):      tgamma=0.85   # decay value for exponential function
150      func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))      def calc_gamma(G,npts):
151      return func          func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))
152            return func
153  gleft  = calc_gamma(tgamma,bleft)
154  gright = calc_gamma(tgamma,bleft)      gleft  = calc_gamma(tgamma,bleft)
155  gbottom= calc_gamma(tgamma,ystep*bn)      gright = calc_gamma(tgamma,bleft)
156        gbottom= calc_gamma(tgamma,ystep*bn)
157  print('gamma', gleft,gright,gbottom)
158        print('gamma', gleft,gright,gbottom)
159  # calculate decay functions
160  def abc_bfunc(gamma,loc,x,G):      # calculate decay functions
161      func=exp(-1.*(gamma*abs(loc-x))**2.)      def abc_bfunc(gamma,loc,x,G):
162      return func          func=exp(-1.*(gamma*abs(loc-x))**2.)
163            return func
164  fleft=abc_bfunc(gleft,bleft,x[0],tgamma)
165  fright=abc_bfunc(gright,bright,x[0],tgamma)      fleft=abc_bfunc(gleft,bleft,x[0],tgamma)
166  fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma)      fright=abc_bfunc(gright,bright,x[0],tgamma)
167  # apply these functions only where relevant      fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma)
168  abcleft=fleft*whereNegative(left)      # apply these functions only where relevant
169  abcright=fright*wherePositive(right)      abcleft=fleft*whereNegative(left)
170  abcbottom=fbottom*wherePositive(bottom)      abcright=fright*wherePositive(right)
171  # make sure the inside of the abc is value 1      abcbottom=fbottom*wherePositive(bottom)
172  abcleft=abcleft+whereZero(abcleft)      # make sure the inside of the abc is value 1
173  abcright=abcright+whereZero(abcright)      abcleft=abcleft+whereZero(abcleft)
174  abcbottom=abcbottom+whereZero(abcbottom)      abcright=abcright+whereZero(abcright)
175  # multiply the conditions together to get a smooth result      abcbottom=abcbottom+whereZero(abcbottom)
176  abc=abcleft*abcright*abcbottom      # multiply the conditions together to get a smooth result
177        abc=abcleft*abcright*abcbottom
178  #visualise the boundary function
179  #abcT=abc.toListOfTuples()      #visualise the boundary function
180  #abcT=np.reshape(abcT,(ndx+1,ndy+1))      #abcT=abc.toListOfTuples()
181  #pl.clf(); pl.imshow(abcT); pl.colorbar();      #abcT=np.reshape(abcT,(ndx+1,ndy+1))
182  #pl.savefig(os.path.join(savepath,"abc.png"))      #pl.clf(); pl.imshow(abcT); pl.colorbar();
183        #pl.savefig(os.path.join(savepath,"abc.png"))
184
185  ############################################FIRST TIME STEPS AND SOURCE
186  # define small radius around point xc      ############################################FIRST TIME STEPS AND SOURCE
187  src_length = 40; print("src_length = ",src_length)      # define small radius around point xc
188  # set initial values for first two time steps with source terms      src_length = 40; print("src_length = ",src_length)
189  y=source[0]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)      # set initial values for first two time steps with source terms
190  src_dir=numpy.array([0.,1.]) # defines direction of point source as down      y=source[0]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)
191  y=y*src_dir      src_dir=numpy.array([0.,1.]) # defines direction of point source as down
192  mypde.setValue(y=y) #set the source as a function on the boundary      y=y*src_dir
193  # turn lumping on for more efficient solving      mypde.setValue(y=y) #set the source as a function on the boundary
194  mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)      # turn lumping on for more efficient solving
195  # for first two time steps      mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
196  u=[0.0,0.0]*wherePositive(x)      # for first two time steps
197  u_m1=u      u=[0.0,0.0]*wherePositive(x)
198        u_m1=u
199  ####################################################ITERATION VARIABLES
200  n=0 # iteration counter      ####################################################ITERATION VARIABLES
201  t=0 # time counter      n=0 # iteration counter
202  ##############################################################ITERATION      t=0 # time counter
203  while t<tend:      ##############################################################ITERATION
204      # get current stress      while t<tend:
205      g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))          # get current stress
206      mypde.setValue(X=-stress*abc) # set PDE values          g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))
207      accel = mypde.getSolution() #get PDE solution for accelleration          mypde.setValue(X=-stress*abc) # set PDE values
208      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement          accel = mypde.getSolution() #get PDE solution for accelleration
209      u_p1=u_p1*abc       # apply boundary conditions          u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
210      u_m1=u; u=u_p1 # shift values by 1          u_p1=u_p1*abc       # apply boundary conditions
211      # save current displacement, acceleration and pressure          u_m1=u; u=u_p1 # shift values by 1
212      if (t >= rtime):          # save current displacement, acceleration and pressure
213          saveVTK(os.path.join(savepath,"ex08b.%05d.vtu"%n),displacement=length(u),\          if (t >= rtime):
214                                      acceleration=length(accel),tensor=stress)              saveVTK(os.path.join(savepath,"ex08b.%05d.vtu"%n),displacement=length(u),\
215          rtime=rtime+rtime_inc #increment data save time                                          acceleration=length(accel),tensor=stress)
216      # increment loop values              rtime=rtime+rtime_inc #increment data save time
217      t=t+h; n=n+1          # increment loop values
218      if (n < ls):          t=t+h; n=n+1
219          y=source[n]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)          if (n < ls):
220          y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the boundary              y=source[n]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)
221      print("time step %d, t=%s"%(n,t))              y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the boundary
222            print("time step %d, t=%s"%(n,t))

Legend:
 Removed from v.5287 changed lines Added in v.5288