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

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

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

revision 3389 by ahallam, Wed Dec 1 23:03:35 2010 UTC revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
1    
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2009-2010 by University of Queensland  # Copyright (c) 2009-2012 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2009-2010 by University of Queensland  __copyright__="""Copyright (c) 2009-2012 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
# Line 45  from esys.finley import ReadMesh Line 46  from esys.finley import ReadMesh
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
# Line 68  vel=Scalar(0,Function(domain)) Line 69  vel=Scalar(0,Function(domain))
69  mu=Scalar(0,Function(domain))  mu=Scalar(0,Function(domain))
70  lam=Scalar(0,Function(domain))  lam=Scalar(0,Function(domain))
71    
72  print 0.5*np.sqrt(prho/(plam+2*pmu))*0.5  print(0.5*np.sqrt(prho/(plam+2*pmu))*0.5)
73    
74  for i in range(0,nlayers):  for i in range(0,nlayers):
75      rho.setTaggedValue('lblock%d'%i,prho[i])      rho.setTaggedValue('lblock%d'%i,prho[i])
# Line 89  lam.setTaggedValue('fault',plam[i]) Line 90  lam.setTaggedValue('fault',plam[i])
90  # Time related variables.  # Time related variables.
91  testing=False  testing=False
92  if testing:  if testing:
93      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.')
94      print "Try changing testing to False for more iterations."      print("Try changing testing to False for more iterations.")
95      tend=0.1      tend=0.1
96  else:  else:
97      tend=0.1    # end time      tend=0.1    # end time
# Line 100  h=0.00001    # time step Line 101  h=0.00001    # time step
101  rtime=0.0 # first time to record  rtime=0.0 # first time to record
102  rtime_inc=tend/750.0 # time increment to record  rtime_inc=tend/750.0 # time increment to record
103  #Check to make sure number of time steps is not too large.  #Check to make sure number of time steps is not too large.
104  print "Time step size= ",h, "Expected number of outputs= ",tend/h  print("Time step size= ",h, "Expected number of outputs= ",tend/h)
105    
106  U0=0.1 # amplitude of point source  U0=0.1 # amplitude of point source
107  ls=500   # length of the source  ls=500   # length of the source
# Line 115  a = 2.0 * (np.pi * dfeq)**2.0 Line 116  a = 2.0 * (np.pi * dfeq)**2.0
116  t0 = 5.0 / (2.0 * np.pi * dfeq)  t0 = 5.0 / (2.0 * np.pi * dfeq)
117  srclength = 5. * t0  srclength = 5. * t0
118  ls = int(srclength/h)  ls = int(srclength/h)
119  print 'source length',ls  print('source length',ls)
120  source=np.zeros(ls,'float') # source array  source=np.zeros(ls,'float') # source array
121  ampmax=0  ampmax=0
122  for it in range(0,ls):  for it in range(0,ls):
# Line 140  mypde.setValue(D=rho*kmat) #set the gene Line 141  mypde.setValue(D=rho*kmat) #set the gene
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 = 10; print "src_length = ",src_length  src_length = 10; print("src_length = ",src_length)
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  yx=(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-xc)-src_length)  yx=(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-xc)-src_length)
# Line 176  while t<tend: Line 177  while t<tend:
177      t=t+h; n=n+1      t=t+h; n=n+1
178      if (n < ls):      if (n < ls):
179          mypde.setValue(y=source[n]*yx*src_dir*stop) #set the source as a function on the boundary          mypde.setValue(y=source[n]*yx*src_dir*stop) #set the source as a function on the boundary
180      print n,"-th time step t ",t      print(n,"-th time step t ",t)

Legend:
Removed from v.3389  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26