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

revision 3346 by caltinay, Fri Nov 12 01:19:02 2010 UTC revision 4576 by sshaw, Mon Dec 9 23:35:30 2013 UTC
# Line 1  Line 1
1
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2009-2010 by University of Queensland  # Copyright (c) 2009-2013 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
# http://www.uq.edu.au/esscc
6  #  #
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
http://www.uq.edu.au/esscc
27  # 3D model with multiple layers.  # 3D model with multiple layers.
28
29  #######################################################EXTERNAL MODULES  #######################################################EXTERNAL MODULES
30    import matplotlib
31    matplotlib.use('agg') #It's just here for automated testing
32  from esys.escript import *  from esys.escript import *
33  from esys.finley import Rectangle  from esys.finley import Rectangle
34  import os  import os
# Line 33  import os Line 36  import os
36  from esys.escript.pdetools import Projector, Locator  from esys.escript.pdetools import Projector, Locator
37  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
38  import numpy as np  import numpy as np
import matplotlib
matplotlib.use('agg') #It's just here for automated testing
39
40  import pylab as pl  import pylab as pl
41  import matplotlib.cm as cm  import matplotlib.cm as cm
# Line 44  from esys.weipa import saveVTK Line 45  from esys.weipa import saveVTK
45
46  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
47  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
48      import sys          import sys
49      print "This example will not run in an MPI world."          print("This example will not run in an MPI world.")
50      sys.exit(0)          sys.exit(0)
51
52  #################################################ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
53  # where to save output data  # where to save output data
54  savepath = "data/example09"  savepath = "data/example09a"
55    meshpath = "data/example09m"
56  mkDir(savepath)  mkDir(savepath)
57  #Geometric and material property related variables.  #Geometric and material property related variables.
mx = 200. # model lenght
my = 200. # model width
mz=100.0
58  step=10.0 # the element size  step=10.0 # the element size
ndx = int(mx/step) # steps in x direction
ndy = int(my/step) # steps in y direction
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/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus  mu2=vel2**2.*rho2/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus
63  lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant  lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant
64
65  # Time related variables.  ####################################################TESTING SWITCH
66  testing=True  testing=True
67  if testing:  if testing:
68      print 'The testing end time is curerntly sellected this severely limits the number of time iterations.'      print('The testing end time is currently selected. This severely limits the number of time iterations.')
69      print "Try changing testing to False for more iterations."      print("Try changing testing to False for more iterations.")
70      tend=0.001      tend=0.001
71        #Model Parameters
72        mx=40.
73        my=40.
74        mz=20.
75        outputs=5
76  else:  else:
77      tend=0.1    # end time      tend=0.1    # end time
78        #Model Parameters
79        mx=100.0   #x width of model
80        my=100.0   #y width of model
81        mz=50.0   #depth of model
82        outputs=200
83
84  h=0.0001    # time step  ####################################################TIME RELATED VARIABLES
85    h=0.00005    # time step
86  # data recording times  # data recording times
87  rtime=0.0 # first time to record  rtime=0.0 # first time to record
88  rtime_inc=tend/50.0 # time increment to record  rtime_inc=tend/outputs # time increment to record
89  #Check to make sure number of time steps is not too large.  #Check to make sure number of time steps is not too large.
90  print "Time step size= ",h, "Expected number of outputs= ",tend/h  print("Time step size= ",h, "Expected number of outputs= ",tend/h)
91
92    ####################################################CREATING THE SOURCE FUNCTION
93  U0=0.1 # amplitude of point source  U0=0.1 # amplitude of point source
94  ls=500   # length of the source  ls=500   # length of the source
95  source=np.zeros(ls,'float') # source array  source=np.zeros(ls,'float') # source array
# Line 95  a = 2.0 * (np.pi * dfeq)**2.0 Line 103  a = 2.0 * (np.pi * dfeq)**2.0
103  t0 = 5.0 / (2.0 * np.pi * dfeq)  t0 = 5.0 / (2.0 * np.pi * dfeq)
104  srclength = 5. * t0  srclength = 5. * t0
105  ls = int(srclength/h)  ls = int(srclength/h)
106  print 'source length',ls  print('source length',ls)
107  source=np.zeros(ls,'float') # source array  source=np.zeros(ls,'float') # source array
108  ampmax=0  ampmax=0
109  for it in range(0,ls):  for it in range(0,ls):
# Line 106  for it in range(0,ls): Line 114  for it in range(0,ls):
114      if (abs(source[it]) > ampmax):      if (abs(source[it]) > ampmax):
115          ampmax = abs(source[it])          ampmax = abs(source[it])
116      time[t]=t*h      time[t]=t*h
#tdecay=decay1*decay2*U0
#decay1=decay1*U0; decay2=decay2*U0
#pl.clf();
#pl.plot(source)
#pl.plot(time,decay1);pl.plot(time,decay2);
#pl.plot(time,tdecay)
#pl.savefig(os.path.join(savepath,'source.png'))
117
118  # will introduce a spherical source at middle left of bottom face  # will introduce a spherical source at middle left of bottom face
119  xc=[50,50,0]  xc=[mx/2,my/2,0]
120
121  ####################################################DOMAIN CONSTRUCTION  ####################################################DOMAIN CONSTRUCTION
123  x=domain.getX() # get the locations of the nodes in the domain  x=domain.getX() # get the locations of the nodes in the domain
124
125  lam=Scalar(0,Function(domain))  lam=Scalar(0,Function(domain))
# Line 135  rho.setTaggedValue("vintfb",rho2) Line 136  rho.setTaggedValue("vintfb",rho2)
136  mypde=LinearPDE(domain) # create pde  mypde=LinearPDE(domain) # create pde
137  mypde.setSymmetryOn() # turn symmetry on  mypde.setSymmetryOn() # turn symmetry on
138  # turn lumping on for more efficient solving  # turn lumping on for more efficient solving
139  #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)  #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING)
140  kmat = kronecker(domain) # create the kronecker delta function of the domain  kmat = kronecker(domain) # create the kronecker delta function of the domain
141  mypde.setValue(D=rho*kmat) #set the general form value D  mypde.setValue(D=rho*kmat) #set the general form value D
142
143  ############################################FIRST TIME STEPS AND SOURCE  ############################################FIRST TIME STEPS AND SOURCE
144  # define small radius around point xc  # define small radius around point xc
146  # set initial values for first two time steps with source terms  # set initial values for first two time steps with source terms
147  xb=FunctionOnBoundary(domain).getX()  xb=FunctionOnBoundary(domain).getX()
149  y=Vector(0.0,FunctionOnBoundary(domain))  stop=Scalar(0.0,FunctionOnBoundary(domain))
150    stop.setTaggedValue("stop",1.0)
151  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
152
153  #sy=sy*src_dir  mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
#sy.setTaggedValue("stop")
y.setTaggedValue("stop",src_dir*source[0])#)*(cos(length(xb-xc)*3.1415/src_length)+1))
mypde.setValue(y=y) #set the source as a function on the boundary
154  # initial value of displacement at point source is constant (U0=0.01)  # initial value of displacement at point source is constant (U0=0.01)
155  # for first two time steps  # for first two time steps
156  u=[0.0,0.0,0.0]*wherePositive(x)  u=[0.0,0.0,0.0]*wherePositive(x)
# Line 163  n=0 # iteration counter Line 161  n=0 # iteration counter
161  t=0 # time counter  t=0 # time counter
162  ##############################################################ITERATION  ##############################################################ITERATION
163  while t<tend:  while t<tend:
164      # get current stress      # get current stress
166      mypde.setValue(X=-stress) # set PDE values      mypde.setValue(X=-stress) # set PDE values
167      accel = mypde.getSolution() #get PDE solution for accelleration      accel = mypde.getSolution() #get PDE solution for accelleration
168      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
169      u_p1=u_p1#*abc          # apply boundary conditions      u_p1=u_p1#*abc          # apply boundary conditions
170      u_m1=u; u=u_p1 # shift values by 1      u_m1=u; u=u_p1 # shift values by 1
171      # save current displacement, acceleration and pressure      # save current displacement, acceleration and pressure
172      if (t >= rtime):      if (t >= rtime):
173          saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\          saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
174                                      acceleration=length(accel),tensor=stress)                                      acceleration=length(accel),tensor=stress)
175          rtime=rtime+rtime_inc #increment data save time          rtime=rtime+rtime_inc #increment data save time
176      # increment loop values      # increment loop values
177      t=t+h; n=n+1      t=t+h; n=n+1
178      if (n < ls):      if (n < ls):
179          y.setTaggedValue("stop",source[n]*src_dir)          mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
180          mypde.setValue(y=y) #set the source as a function on the boundary      print("time step %d, t=%s"%(n,t))
print n,"-th time step t ",t

Legend:
 Removed from v.3346 changed lines Added in v.4576