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

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

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

revision 3387 by caltinay, Thu Nov 25 07:09:23 2010 UTC revision 3892 by jfenwick, Tue Apr 10 08:57:23 2012 UTC
# Line 51  import numpy as np Line 51  import numpy as np
51  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
52  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
53      import sys      import sys
54      print "This example will not run in an MPI world."      print("This example will not run in an MPI world.")
55      sys.exit(0)      sys.exit(0)
56    
57  #################################################ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
# Line 87  lam2=mu2*6.; lam1=mu1*6. #lames constant Line 87  lam2=mu2*6.; lam1=mu1*6. #lames constant
87  # Time related variables.  # Time related variables.
88  testing=True  testing=True
89  if testing:  if testing:
90      print 'The testing end time is currently selected. This severely limits the number of time iterations.'      print('The testing end time is currently selected. This severely limits the number of time iterations.')
91      print "Try changing testing to False for more iterations."      print("Try changing testing to False for more iterations.")
92      tend=0.001      tend=0.001
93  else:  else:
94      tend=0.5    # end time      tend=0.5    # end time
# Line 100  rtime_inc=tend/50.0 # time increment to Line 100  rtime_inc=tend/50.0 # time increment to
100  # will introduce a spherical source at middle left of bottom face  # will introduce a spherical source at middle left of bottom face
101  xc=[width/2,0]  xc=[width/2,0]
102  #Check to make sure number of time steps is not too large.  #Check to make sure number of time steps is not too large.
103  print "Time step size= ",h, "Expected number of outputs= ",tend/h  print("Time step size= ",h, "Expected number of outputs= ",tend/h)
104    
105  U0=0.1 # amplitude of point source  U0=0.1 # amplitude of point source
106  ls=500   # length of the source  ls=500   # length of the source
# Line 115  a = 2.0 * (np.pi * dfeq)**2.0 Line 115  a = 2.0 * (np.pi * dfeq)**2.0
115  t0 = 5.0 / (2.0 * np.pi * dfeq)  t0 = 5.0 / (2.0 * np.pi * dfeq)
116  srclength = 5. * t0  srclength = 5. * t0
117  ls = int(srclength/h)  ls = int(srclength/h)
118  print 'source length',ls  print('source length',ls)
119  source=np.zeros(ls,'float') # source array  source=np.zeros(ls,'float') # source array
120  ampmax=0  ampmax=0
121  for it in range(0,ls):  for it in range(0,ls):
# Line 181  d.setScriptFileName(os.path.join(save_pa Line 181  d.setScriptFileName(os.path.join(save_pa
181  d.setMeshFileName(os.path.join(save_path,"example08c.msh"))  d.setMeshFileName(os.path.join(save_path,"example08c.msh"))
182  domain=MakeDomain(d, optimizeLabeling=True)  domain=MakeDomain(d, optimizeLabeling=True)
183  x=domain.getX()  x=domain.getX()
184  print "Domain has been generated ..."  print("Domain has been generated ...")
185    
186  lam=Scalar(0,Function(domain))  lam=Scalar(0,Function(domain))
187  lam.setTaggedValue("top",lam1)  lam.setTaggedValue("top",lam1)
# Line 219  gleft  = calc_gamma(tgamma,bleft) Line 219  gleft  = calc_gamma(tgamma,bleft)
219  gright = calc_gamma(tgamma,bleft)  gright = calc_gamma(tgamma,bleft)
220  gbottom= calc_gamma(tgamma,ystep*bn)  gbottom= calc_gamma(tgamma,ystep*bn)
221    
222  print 'gamma', gleft,gright,gbottom  print('gamma', gleft,gright,gbottom)
223    
224  # calculate decay functions  # calculate decay functions
225  def abc_bfunc(gamma,loc,x,G):  def abc_bfunc(gamma,loc,x,G):
# Line 242  abc=abcleft*abcright*abcbottom Line 242  abc=abcleft*abcright*abcbottom
242    
243  ############################################FIRST TIME STEPS AND SOURCE  ############################################FIRST TIME STEPS AND SOURCE
244  # define small radius around point xc  # define small radius around point xc
245  src_length = 40; print "src_length = ",src_length  src_length = 40; print("src_length = ",src_length)
246  # set initial values for first two time steps with source terms  # set initial values for first two time steps with source terms
247  xb=FunctionOnBoundary(domain).getX()  xb=FunctionOnBoundary(domain).getX()
248  y=source[0]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))  y=source[0]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length))
# Line 278  while t<tend: Line 278  while t<tend:
278      if (n < ls):      if (n < ls):
279          y=source[n]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)          y=source[n]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)
280          y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the boundary          y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the boundary
281      print n,"-th time step t ",t      print(n,"-th time step t ",t)

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

  ViewVC Help
Powered by ViewVC 1.1.26