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

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

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

revision 5287 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 1  Line 1 
1  from __future__ import division  from __future__ import division, print_function
 from __future__ import print_function  
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2009-2014 by University of Queensland  # Copyright (c) 2009-2014 by University of Queensland
# Line 32  __url__="https://launchpad.net/escript-f Line 31  __url__="https://launchpad.net/escript-f
31  import matplotlib  import matplotlib
32  matplotlib.use('agg') #It's just here for automated testing  matplotlib.use('agg') #It's just here for automated testing
33  from esys.escript import *  from esys.escript import *
 from esys.finley import Rectangle  
34  from esys.weipa import saveVTK  from esys.weipa import saveVTK
35  import os  import os
36  # smoothing operator  # smoothing operator
# Line 43  import numpy as np Line 41  import numpy as np
41  import pylab as pl  import pylab as pl
42  import matplotlib.cm as cm  import matplotlib.cm as cm
43  from esys.escript.linearPDEs import LinearPDE, SolverOptions  from esys.escript.linearPDEs import LinearPDE, SolverOptions
44  from esys.finley import ReadMesh  try:
45        # This imports the rectangle domain function
46        from esys.finley import Rectangle, ReadMesh
47        HAVE_FINLEY = True
48    except ImportError:
49        print("Finley module not available")
50        HAVE_FINLEY = False
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  if HAVE_FINLEY:
58  # where to save output data      #################################################ESTABLISHING VARIABLES
59  savepath = "data/example09b"      # where to save output data
60  meshpath = "data/example09m"      savepath = "data/example09b"
61  mkDir(savepath)      meshpath = "data/example09m"
62  #Geometric and material property related variables.      mkDir(savepath)
63  step=4.0 # the element size      #Geometric and material property related variables.
64        step=4.0 # the element size
65  vel=1800.    #starting velocity  
66  rhoc=2000.   #starting density      vel=1800.    #starting velocity
67  nlayers=9  #number of layers in layercake model.      rhoc=2000.   #starting density
68        nlayers=9  #number of layers in layercake model.
69  ####################################################TESTING SWITCH  
70  testing=True      ####################################################TESTING SWITCH
71  if testing:      testing=True
72      print('The testing end time is currently selected. This severely limits the number of time iterations.')      if testing:
73      print("Try changing testing to False for more iterations.")          print('The testing end time is currently selected. This severely limits the number of time iterations.')
74      tend=0.001          print("Try changing testing to False for more iterations.")
75      #Model Parameters              tend=0.001
76      mx=40.          #Model Parameters    
77      my=40.          mx=40.
78      mz=20.          my=40.
79      outputs=5          mz=20.
80  else:          outputs=5
81      tend=0.1    # end time      else:
82      #Model Parameters          tend=0.1    # end time
83      mx=100.0   #x width of model          #Model Parameters
84      my=100.0   #y width of model          mx=100.0   #x width of model
85      mz=50.0   #depth of model          my=100.0   #y width of model
86      outputs=200          mz=50.0   #depth of model
87            outputs=200
88  ####################################################TIME RELATED VARIABLES  
89  h=0.00001    # time step      ####################################################TIME RELATED VARIABLES
90  # data recording times      h=0.00001    # time step
91  rtime=0.0 # first time to record      # data recording times
92  rtime_inc=tend/outputs # time increment to record      rtime=0.0 # first time to record
93  #Check to make sure number of time steps is not too large.      rtime_inc=tend/outputs # time increment to record
94  print("Time step size= ",h, "Expected number of outputs= ",tend/h)      #Check to make sure number of time steps is not too large.
95        print("Time step size= ",h, "Expected number of outputs= ",tend/h)
96  ####################################################CREATING THE SOURCE FUNCTION  
97  U0=0.1 # amplitude of point source      ####################################################CREATING THE SOURCE FUNCTION
98  ls=500   # length of the source      U0=0.1 # amplitude of point source
99  source=np.zeros(ls,'float') # source array      ls=500   # length of the source
100  decay1=np.zeros(ls,'float') # decay curve one      source=np.zeros(ls,'float') # source array
101  decay2=np.zeros(ls,'float') # decay curve two      decay1=np.zeros(ls,'float') # decay curve one
102  time=np.zeros(ls,'float')   # time values      decay2=np.zeros(ls,'float') # decay curve two
103  g=np.log(0.01)/ls      time=np.zeros(ls,'float')   # time values
104        g=np.log(0.01)/ls
105  dfeq=50 #Dominant Frequency  
106  a = 2.0 * (np.pi * dfeq)**2.0      dfeq=50 #Dominant Frequency
107  t0 = 5.0 / (2.0 * np.pi * dfeq)      a = 2.0 * (np.pi * dfeq)**2.0
108  srclength = 5. * t0      t0 = 5.0 / (2.0 * np.pi * dfeq)
109  ls = int(srclength/h)      srclength = 5. * t0
110  print('source length',ls)      ls = int(srclength/h)
111  source=np.zeros(ls,'float') # source array      print('source length',ls)
112  ampmax=0      source=np.zeros(ls,'float') # source array
113  for it in range(0,ls):      ampmax=0
114      t = it*h      for it in range(0,ls):
115      tt = t-t0          t = it*h
116      dum1 = np.exp(-a * tt * tt)          tt = t-t0
117      source[it] = -2. * a * tt * dum1          dum1 = np.exp(-a * tt * tt)
118      if (abs(source[it]) > ampmax):          source[it] = -2. * a * tt * dum1
119          ampmax = abs(source[it])          if (abs(source[it]) > ampmax):
120      time[t]=t*h              ampmax = abs(source[it])
121            time[t]=t*h
122  # will introduce a spherical source at middle left of bottom face  
123  xc=[mx/2,my/2,0]      # will introduce a spherical source at middle left of bottom face
124        xc=[mx/2,my/2,0]
125  ####################################################DOMAIN CONSTRUCTION  
126  domain=ReadMesh(os.path.join(meshpath,'example09lc.fly')) # create the domain      ####################################################DOMAIN CONSTRUCTION
127  x=domain.getX() # get the locations of the nodes in the domain      domain=ReadMesh(os.path.join(meshpath,'example09lc.fly')) # create the domain
128        x=domain.getX() # get the locations of the nodes in the domain
129  lam=Scalar(0,Function(domain))  
130  mu=Scalar(0,Function(domain))      lam=Scalar(0,Function(domain))
131  rho=Scalar(0,Function(domain))      mu=Scalar(0,Function(domain))
132        rho=Scalar(0,Function(domain))
133  #Setting parameters for each layer in the model.  
134  for i in range(0,nlayers):      #Setting parameters for each layer in the model.
135      rho.setTaggedValue("volume_%d"%i,rhoc+i*100.)      for i in range(0,nlayers):
136      lamc=(vel+i*100.)**2.*(rhoc+i*100.)/2.          rho.setTaggedValue("volume_%d"%i,rhoc+i*100.)
137      muc=(vel+i*100.)**2.*(rhoc+i*100.)/4.          lamc=(vel+i*100.)**2.*(rhoc+i*100.)/2.
138      lam.setTaggedValue("volume_%d"%i,lamc)          muc=(vel+i*100.)**2.*(rhoc+i*100.)/4.
139      mu.setTaggedValue("volume_%d"%i,muc)          lam.setTaggedValue("volume_%d"%i,lamc)
140            mu.setTaggedValue("volume_%d"%i,muc)
141  ##########################################################ESTABLISH PDE  
142  mypde=LinearPDE(domain) # create pde      ##########################################################ESTABLISH PDE
143  mypde.setSymmetryOn() # turn symmetry on      mypde=LinearPDE(domain) # create pde
144  # turn lumping on for more efficient solving      mypde.setSymmetryOn() # turn symmetry on
145  #mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)      # turn lumping on for more efficient solving
146  kmat = kronecker(domain) # create the kronecker delta function of the domain      #mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
147  mypde.setValue(D=rho*kmat) #set the general form value D      kmat = kronecker(domain) # create the kronecker delta function of the domain
148        mypde.setValue(D=rho*kmat) #set the general form value D
149  ############################################FIRST TIME STEPS AND SOURCE  
150  # define small radius around point xc      ############################################FIRST TIME STEPS AND SOURCE
151  src_rad = 20; print("src radius= ",src_rad)      # define small radius around point xc
152  # set initial values for first two time steps with source terms      src_rad = 20; print("src radius= ",src_rad)
153  xb=FunctionOnBoundary(domain).getX()      # set initial values for first two time steps with source terms
154  yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)      xb=FunctionOnBoundary(domain).getX()
155  stop=Scalar(0.0,FunctionOnBoundary(domain))      yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)
156  stop.setTaggedValue("intface_0",1.0)      stop=Scalar(0.0,FunctionOnBoundary(domain))
157  src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down      stop.setTaggedValue("intface_0",1.0)
158        src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
159  mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary  
160  # initial value of displacement at point source is constant (U0=0.01)      mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
161  # for first two time steps      # initial value of displacement at point source is constant (U0=0.01)
162  u=[0.0,0.0,0.0]*wherePositive(x)      # for first two time steps
163  u_m1=u      u=[0.0,0.0,0.0]*wherePositive(x)
164        u_m1=u
165  ####################################################ITERATION VARIABLES  
166  n=0 # iteration counter      ####################################################ITERATION VARIABLES
167  t=0 # time counter      n=0 # iteration counter
168  ##############################################################ITERATION      t=0 # time counter
169  while t<tend:      ##############################################################ITERATION
170      # get current stress      while t<tend:
171      g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc          # get current stress
172      mypde.setValue(X=-stress) # set PDE values          g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
173      accel = mypde.getSolution() #get PDE solution for accelleration          mypde.setValue(X=-stress) # set PDE values
174      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement          accel = mypde.getSolution() #get PDE solution for accelleration
175      u_p1=u_p1#*abc          # apply boundary conditions          u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
176      u_m1=u; u=u_p1 # shift values by 1          u_p1=u_p1#*abc          # apply boundary conditions
177      # save current displacement, acceleration and pressure          u_m1=u; u=u_p1 # shift values by 1
178      if (t >= rtime):          # save current displacement, acceleration and pressure
179          saveVTK(os.path.join(savepath,"ex09b.%05d.vtu"%n),displacement=length(u),\          if (t >= rtime):
180                                      acceleration=length(accel),tensor=stress)              saveVTK(os.path.join(savepath,"ex09b.%05d.vtu"%n),displacement=length(u),\
181          rtime=rtime+rtime_inc #increment data save time                                          acceleration=length(accel),tensor=stress)
182      # increment loop values              rtime=rtime+rtime_inc #increment data save time
183      t=t+h; n=n+1          # increment loop values
184      if (n < ls):          t=t+h; n=n+1
185          mypde.setValue(y=source[n]*yx*src_dir*stop) #set the source as a function on the boundary          if (n < ls):
186      print("time step %d, t=%s"%(n,t))              mypde.setValue(y=source[n]*yx*src_dir*stop) #set the source as a function on the boundary
187            print("time step %d, t=%s"%(n,t))

Legend:
Removed from v.5287  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26