/[escript]/trunk/doc/examples/cookbook/example09a.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3346 by caltinay, Fri Nov 12 01:19:02 2010 UTC revision 3892 by jfenwick, Tue Apr 10 08:57:23 2012 UTC
# Line 45  from esys.weipa import saveVTK Line 45  from esys.weipa import saveVTK
45  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
46  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
47      import sys      import sys
48      print "This example will not run in an MPI world."      print("This example will not run in an MPI world.")
49      sys.exit(0)      sys.exit(0)
50    
51  #################################################ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
52  # where to save output data  # where to save output data
53  savepath = "data/example09"  savepath = "data/example09a"
54    meshpath = "data/example09m"
55  mkDir(savepath)  mkDir(savepath)
56  #Geometric and material property related variables.  #Geometric and material property related variables.
 mx = 200. # model lenght  
 my = 200. # model width  
 mz=100.0  
57  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)  
58    
59  vel2=1800.;   vel1=3000.  vel2=1800.;   vel1=3000.
60  rho2=2300.;   rho1=3100. #density  rho2=2300.;   rho1=3100. #density
61  mu2=vel2**2.*rho2/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus  mu2=vel2**2.*rho2/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus
62  lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant  lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant
63    
64  # Time related variables.  ####################################################TESTING SWITCH
65  testing=True  testing=True
66  if testing:  if testing:
67      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.')
68      print "Try changing testing to False for more iterations."      print("Try changing testing to False for more iterations.")
69      tend=0.001      tend=0.001
70        #Model Parameters    
71        mx=40.
72        my=40.
73        mz=20.
74        outputs=5
75  else:  else:
76      tend=0.1    # end time      tend=0.1    # end time
77        #Model Parameters
78        mx=100.0   #x width of model
79        my=100.0   #y width of model
80        mz=50.0   #depth of model
81        outputs=200
82    
83  h=0.0001    # time step  ####################################################TIME RELATED VARIABLES
84    h=0.00005    # time step
85  # data recording times  # data recording times
86  rtime=0.0 # first time to record  rtime=0.0 # first time to record
87  rtime_inc=tend/50.0 # time increment to record  rtime_inc=tend/outputs # time increment to record
88  #Check to make sure number of time steps is not too large.  #Check to make sure number of time steps is not too large.
89  print "Time step size= ",h, "Expected number of outputs= ",tend/h  print("Time step size= ",h, "Expected number of outputs= ",tend/h)
90    
91    ####################################################CREATING THE SOURCE FUNCTION
92  U0=0.1 # amplitude of point source  U0=0.1 # amplitude of point source
93  ls=500   # length of the source  ls=500   # length of the source
94  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 102  a = 2.0 * (np.pi * dfeq)**2.0
102  t0 = 5.0 / (2.0 * np.pi * dfeq)  t0 = 5.0 / (2.0 * np.pi * dfeq)
103  srclength = 5. * t0  srclength = 5. * t0
104  ls = int(srclength/h)  ls = int(srclength/h)
105  print 'source length',ls  print('source length',ls)
106  source=np.zeros(ls,'float') # source array  source=np.zeros(ls,'float') # source array
107  ampmax=0  ampmax=0
108  for it in range(0,ls):  for it in range(0,ls):
# Line 106  for it in range(0,ls): Line 113  for it in range(0,ls):
113      if (abs(source[it]) > ampmax):      if (abs(source[it]) > ampmax):
114          ampmax = abs(source[it])          ampmax = abs(source[it])
115      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'))  
116    
117  # will introduce a spherical source at middle left of bottom face  # will introduce a spherical source at middle left of bottom face
118  xc=[50,50,0]  xc=[mx/2,my/2,0]
119    
120  ####################################################DOMAIN CONSTRUCTION  ####################################################DOMAIN CONSTRUCTION
121  domain=ReadMesh(os.path.join(savepath,'example09m.fly')) # create the domain  domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
122  x=domain.getX() # get the locations of the nodes in the domain  x=domain.getX() # get the locations of the nodes in the domain
123    
124  lam=Scalar(0,Function(domain))  lam=Scalar(0,Function(domain))
# Line 135  rho.setTaggedValue("vintfb",rho2) Line 135  rho.setTaggedValue("vintfb",rho2)
135  mypde=LinearPDE(domain) # create pde  mypde=LinearPDE(domain) # create pde
136  mypde.setSymmetryOn() # turn symmetry on  mypde.setSymmetryOn() # turn symmetry on
137  # turn lumping on for more efficient solving  # turn lumping on for more efficient solving
138  #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)  #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING)
139  kmat = kronecker(domain) # create the kronecker delta function of the domain  kmat = kronecker(domain) # create the kronecker delta function of the domain
140  mypde.setValue(D=rho*kmat) #set the general form value D  mypde.setValue(D=rho*kmat) #set the general form value D
141    
142  ############################################FIRST TIME STEPS AND SOURCE  ############################################FIRST TIME STEPS AND SOURCE
143  # define small radius around point xc  # define small radius around point xc
144  src_length = 20; print "src_length = ",src_length  src_rad = 20; print("sourc radius= ",src_rad)
145  # set initial values for first two time steps with source terms  # set initial values for first two time steps with source terms
146  xb=FunctionOnBoundary(domain).getX()  xb=FunctionOnBoundary(domain).getX()
147  #sy=source[0]*(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))  yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)
148  y=Vector(0.0,FunctionOnBoundary(domain))  stop=Scalar(0.0,FunctionOnBoundary(domain))
149    stop.setTaggedValue("stop",1.0)
150  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
151    
152  #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  
153  # initial value of displacement at point source is constant (U0=0.01)  # initial value of displacement at point source is constant (U0=0.01)
154  # for first two time steps  # for first two time steps
155  u=[0.0,0.0,0.0]*wherePositive(x)  u=[0.0,0.0,0.0]*wherePositive(x)
# Line 172  while t<tend: Line 169  while t<tend:
169      u_m1=u; u=u_p1 # shift values by 1      u_m1=u; u=u_p1 # shift values by 1
170      # save current displacement, acceleration and pressure      # save current displacement, acceleration and pressure
171      if (t >= rtime):      if (t >= rtime):
172          saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\          saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
173                                      acceleration=length(accel),tensor=stress)                                      acceleration=length(accel),tensor=stress)
174          rtime=rtime+rtime_inc #increment data save time          rtime=rtime+rtime_inc #increment data save time
175      # increment loop values      # increment loop values
176      t=t+h; n=n+1      t=t+h; n=n+1
177      if (n < ls):      if (n < ls):
178          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
179          mypde.setValue(y=y) #set the source as a function on the boundary      print(n,"-th time step t ",t)
     print n,"-th time step t ",t  

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

  ViewVC Help
Powered by ViewVC 1.1.26