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

revision 3025 by ahallam, Thu May 6 01:20:46 2010 UTC revision 3029 by ahallam, Fri May 21 02:01:37 2010 UTC
21
23  # example07a.py  # example08a.py
24  # Antony Hallam  # Antony Hallam
25  # Seismic Wave Equation Simulation using acceleration solution.  # Seismic Wave Equation Simulation using acceleration solution.
# You can shorten the execution time by reducing variable tend from 60 to 0.5
26
27  #######################################################EXTERNAL MODULES  #######################################################EXTERNAL MODULES
28  from esys.escript import *  from esys.escript import *
# Line 40  from esys.escript.linearPDEs import Line Line 39  from esys.escript.linearPDEs import Line
39
40  #################################################ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
41  # where to save output data  # where to save output data
42  savepath = "data/example07a"  savepath = "data/example08a"
43  mkDir(savepath)  mkDir(savepath)
44  #Geometric and material property related variables.  #Geometric and material property related variables.
45  mx = 1000. # model lenght  mx = 1000. # model lenght
# Line 55  rho=1154.   #density Line 54  rho=1154.   #density
54  # Time related variables.  # Time related variables.
55  tend=1.5    # end time  tend=1.5    # end time
56  h=0.001     # time step  h=0.001     # time step
57    # data recording times
58    rtime=0.0 # first time to record
59    rtime_inc=tend/20.0 # time increment to record
60  #Check to make sure number of time steps is not too large.  #Check to make sure number of time steps is not too large.
61  print "Time step size= ",h, "Expected number of outputs= ",tend/h  print "Time step size= ",h, "Expected number of outputs= ",tend/h
62
# Line 69  x=domain.getX() # get the locations of t Line 71  x=domain.getX() # get the locations of t
71  ##########################################################ESTABLISH PDE  ##########################################################ESTABLISH PDE
72  mypde=LinearPDE(domain) # create pde  mypde=LinearPDE(domain) # create pde
73  mypde.setSymmetryOn() # turn symmetry on  mypde.setSymmetryOn() # turn symmetry on
mypde.setValue(D=1.) # set the value of D in the general form to 1.
74  # turn lumping on for more efficient solving  # turn lumping on for more efficient solving
75  mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)  mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
76  kmat = kronecker(domain) # create the kronecker delta function of the domain  kmat = kronecker(domain) # create the kronecker delta function of the domain
# Line 79  mypde.setValue(D=kmat*rho) #set the gene Line 80  mypde.setValue(D=kmat*rho) #set the gene
80  # define small radius around point xc  # define small radius around point xc
81  src_length = 20; print "src_length = ",src_length  src_length = 20; print "src_length = ",src_length
82  # set initial values for first two time steps with source terms  # set initial values for first two time steps with source terms
83    y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)
84  src_dir=numpy.array([0.,-1.]) # defines direction of point source as down  src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
85  y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)*src_dir  y=y*src_dir
86  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
87  # initial value of displacement at point source is constant (U0=0.01)  # initial value of displacement at point source is constant (U0=0.01)
88  # for first two time steps  # for first two time steps
# Line 100  while t<tend: Line 102  while t<tend:
102      u_m1=u; u=u_p1 # shift values by 1      u_m1=u; u=u_p1 # shift values by 1
103      # save current displacement, acceleration and pressure      # save current displacement, acceleration and pressure
104      if (t >= rtime):      if (t >= rtime):
105          saveVTK(os.path.join(savepath,"ex07b.%i.vtu"%n),displacement=length(u),\          saveVTK(os.path.join(savepath,"ex08a.%05d.vtu"%n),displacement=length(u),\
106                                      acceleration=length(accel),tensor=pres)                                      acceleration=length(accel),tensor=stress)
107          rtime=rtime+rtime_inc #increment data save time          rtime=rtime+rtime_inc #increment data save time
108      # increment loop values      # increment loop values
109      t=t+h; n=n+1      t=t+h; n=n+1

Legend:
 Removed from v.3025 changed lines Added in v.3029