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

  ViewVC Help
Powered by ViewVC 1.1.26