/[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 4087 by caltinay, Thu Nov 22 22:28:01 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 26  __url__="https://launchpad.net/escript-f Line 27  __url__="https://launchpad.net/escript-f
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 45  from esys.weipa import saveVTK Line 46  from esys.weipa import saveVTK
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
122  domain=ReadMesh(os.path.join(savepath,'example09m.fly')) # create the domain  domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
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
145  src_length = 20; print "src_length = ",src_length  src_rad = 20; print("sourc radius= ",src_rad)
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()
148  #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)
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 172  while t<tend: Line 170  while t<tend:
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.4087

  ViewVC Help
Powered by ViewVC 1.1.26